1. Obtenemos los paquetes a usar

#paquetes <- c("cluster", "dbscan", "e1071", "factoextra", "fpc", "GGally", "kernlab", "mclust", "mlbench", "scatterpie", "seriation", "tidyverse")

#install.packages(paquetes, dependencies = TRUE)

¿Para qué sirven estos paquetes?

  1. “cluster”: Este paquete proporciona una variedad de funciones para la agrupación de datos, incluyendo el agrupamiento jerárquico y k-means.
  2. “dbscan”: Este paquete proporciona funciones para realizar el agrupamiento basado en densidad DBSCAN.
  3. “e1071”: Este paquete proporciona una variedad de funciones para aprendizaje automático, incluyendo funciones para clasificación, análisis discriminante y clasificación bayesiana.
  4. “factoextra”: Este paquete proporciona funciones adicionales para realizar análisis factorial, incluyendo gráficos y visualizaciones.
  5. “fpc”: Este paquete proporciona funciones para el agrupamiento basado en densidad y la validación de modelos de agrupamiento.
  6. “GGally”: Este paquete proporciona funciones para crear gráficos de pareja y matrices de scatter plot, con la librería ggplot2.
  7. “kernlab”: Este paquete proporciona funciones para aprendizaje automático basado en kernel, incluyendo funciones para clasificación, regresión y agrupamiento.
  8. “mclust”: Este paquete proporciona funciones para agrupamiento basado en modelos y validación de modelos de agrupamiento.
  9. “mlbench”: Este paquete proporciona conjuntos de datos y funciones para el aprendizaje automático, incluyendo conjuntos de datos para clasificación, regresión y agrupamiento.
  10. “scatterpie”: Este paquete proporciona funciones para crear gráficos de dispersión con pie charts en cada punto.
  11. “seriation”: Este paquete proporciona funciones para realizar la seriación, una técnica de ordenamiento de datos.
  12. “tidyverse”: Este paquete es un conjunto de paquetes de R diseñados para trabajar juntos para facilitar la manipulación, visualización y análisis de datos en R. Incluye paquetes como ggplot2, dplyr, tidyr, entre otros.

2. Importamos los datos (data preparation)

Abrimos Tidyverse

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

Abrimos la base de datos de la población penitenciaria

library(readxl)

# Lee la hoja por defecto del archivo Excel
datos <- read_excel("Denuncias.xlsx")
summary (datos)
##   Ubigeo Hecho    Departamento        Provincia           Distrito        
##  Min.   : 10101   Length:1435        Length:1435        Length:1435       
##  1st Qu.: 50621   Class :character   Class :character   Class :character  
##  Median :110210   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :112629                                                           
##  3rd Qu.:151030                                                           
##  Max.   :250401                                                           
##                                                                           
##      ESTAFA          EXTORSION          HURTO             ROBO        
##  Min.   :   1.00   Min.   :   1.0   Min.   :   1.0   Min.   :   1.00  
##  1st Qu.:   2.00   1st Qu.:   1.0   1st Qu.:   2.0   1st Qu.:   2.00  
##  Median :   5.00   Median :   4.5   Median :   8.0   Median :   4.00  
##  Mean   :  40.31   Mean   :  31.2   Mean   : 143.4   Mean   :  95.77  
##  3rd Qu.:  20.00   3rd Qu.:  19.0   3rd Qu.:  42.0   3rd Qu.:  24.00  
##  Max.   :1278.00   Max.   :1092.0   Max.   :9260.0   Max.   :4428.00  
##  NA's   :646       NA's   :869      NA's   :123      NA's   :389      
##    HOMICIDIO      
##  Min.   :  1.000  
##  1st Qu.:  1.000  
##  Median :  2.000  
##  Mean   :  5.265  
##  3rd Qu.:  4.000  
##  Max.   :131.000  
##  NA's   :1050
head (datos)
## # A tibble: 6 × 9
##   `Ubigeo Hecho` Departamento Provincia   Distrito  ESTAFA EXTORSION HURTO  ROBO
##            <dbl> <chr>        <chr>       <chr>      <dbl>     <dbl> <dbl> <dbl>
## 1          10101 AMAZONAS     CHACHAPOYAS CHACHAPO…    127        34   172    31
## 2          10103 AMAZONAS     CHACHAPOYAS BALSAS        NA         1     3     1
## 3          10105 AMAZONAS     CHACHAPOYAS CHILIQUIN     NA        NA    NA     1
## 4          10106 AMAZONAS     CHACHAPOYAS CHUQUIBA…     NA        NA     2    NA
## 5          10109 AMAZONAS     CHACHAPOYAS LA JALCA       3        NA     1     1
## 6          10110 AMAZONAS     CHACHAPOYAS LEIMEBAM…      1        NA     6    NA
## # ℹ 1 more variable: HOMICIDIO <dbl>
tail (datos)
## # A tibble: 6 × 9
##   `Ubigeo Hecho` Departamento Provincia  Distrito   ESTAFA EXTORSION HURTO  ROBO
##            <dbl> <chr>        <chr>      <chr>       <dbl>     <dbl> <dbl> <dbl>
## 1         250301 UCAYALI      PADRE ABAD PADRE ABAD     16         7    85    14
## 2         250302 UCAYALI      PADRE ABAD IRAZOLA         6         2    27    11
## 3         250303 UCAYALI      PADRE ABAD CURIMANA        1        NA    27     3
## 4         250304 UCAYALI      PADRE ABAD NESHUYA         3         4    36     7
## 5         250305 UCAYALI      PADRE ABAD ALEXANDER…      2        NA    12     4
## 6         250401 UCAYALI      PURÚS      PURÚS           2        NA     9    NA
## # ℹ 1 more variable: HOMICIDIO <dbl>
library(dplyr)

# Eliminar filas con NA en cualquier columna
data_limpia <- datos %>% drop_na()

3. Métodos de conglomeración (clustering methods)

3.1 Conglomerados k-medias (k-means Clustering)

# Seleccionamos las variables relevantes (EDAD y Nº DE INGRESOS)
variables_numericas <- data_limpia[, c("ESTAFA", "EXTORSION", "HURTO", "ROBO", "HOMICIDIO")]

# Escalamos las variables
variables_escaladas <- scale(variables_numericas)
# Aplicamos k-means con 4 conglomerados y 10 inicios aleatorios
set.seed(123)
km <- kmeans(variables_numericas, centers = 3, nstart = 10)

# Mostramos los resultados
km
## K-means clustering with 3 clusters of sizes 13, 215, 43
## 
## Cluster means:
##      ESTAFA EXTORSION     HURTO       ROBO HOMICIDIO
## 1 724.76923 405.00000 4830.3077 2800.23077  47.61538
## 2  38.48837  29.20465  187.8791   90.91163   3.60000
## 3 241.46512 110.81395 1513.4884  868.93023  11.32558
## 
## Clustering vector:
##   [1] 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 3 2 2 3 2 2 2 2 3 2 2 2 2 3 2 2 2 2 2 3 2
##  [38] 2 2 2 2 3 2 2 2 2 2 3 2 2 1 2 2 2 2 3 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2
## [112] 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 3 3 2 2 2 2 2 2 2 2 2 2 2
## [149] 2 2 2 1 2 1 3 3 3 2 3 2 1 3 3 3 3 1 3 1 3 2 2 1 2 2 3 3 3 3 1 3 1 3 3 2 1
## [186] 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [223] 1 3 2 2 2 2 3 2 2 2 2 2 2 2 3 2 2 2 3 2 3 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2
## [260] 2 2 2 2 2 2 2 3 2 2 2 2
## 
## Within cluster sum of squares by cluster:
## [1] 58241224 11614616 23279782
##  (between_SS / total_SS =  81.8 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
# Creamos un data frame escalado con la asignación de clusters
clustered_data <- as.data.frame(variables_numericas) %>%
  add_column(cluster = factor(km$cluster))

# Visualizamos las primeras filas
head(clustered_data)
##   ESTAFA EXTORSION HURTO ROBO HOMICIDIO cluster
## 1    127        34   172   31         2       2
## 2      2         1     8    5         4       2
## 3     86        16   552  203         9       2
## 4     17         1    69   14         2       2
## 5      1         1     8    1         1       2
## 6     43        18   189   56         1       2
library(ggplot2)

# Visualizamos los clusters
ggplot(clustered_data, aes(x = HURTO, y = ROBO, color = cluster)) +
  geom_point(size = 2) +
  labs(title = "Distribución de Conglomerados k-means", x = "Hurto", y = "Robo") +
  theme_minimal()

# Visualizamos con 8 clusters
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
fviz_cluster(kmeans(variables_numericas, centers = 3, nstart = 10),
             data = variables_numericas, centroids = TRUE, geom = "point", ellipse.type = "norm") +
  labs(title = "Distribución con 4 Clusters")

#3.2. Análisis Jerárquico

# Seleccionamos una muestra de 10,000 registros (ajustar según la memoria disponible)
set.seed(123)  # Semilla para reproducibilidad

# Creamos la matriz de distancias y continuar con el análisis
d <- dist(variables_escaladas)
# Agrupación jerárquica aglomerativa con enlace completo
hc <- hclust(d, method = "complete")

# Visualizamos el dendrograma
plot(hc, main = "Dendrograma - Método Complete", xlab = "", sub = "", cex = 0.8)

library(factoextra)

# Visualizamos el dendrograma cortado en 4 conglomerados
fviz_dend(hc, k = 3, rect = TRUE, show_labels = FALSE, cex = 0.6)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <]8;;https://github.com/kassambara/factoextra/issueshttps://github.com/kassambara/factoextra/issues]8;;>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Extraemos los clusters
clusters <- cutree(hc, k = 3)

# Creamos un data frame con las asignaciones
cluster_complete <- as.data.frame(variables_escaladas) %>%
  add_column(cluster = factor(clusters))

# Vemos los primeros registros con clusters asignados
head(cluster_complete)
##        ESTAFA  EXTORSION       HURTO       ROBO  HOMICIDIO cluster
## 1  0.13158691 -0.2052215 -0.38528446 -0.4494547 -0.3489628       1
## 2 -0.57182267 -0.4638960 -0.52603857 -0.4867504 -0.2076042       1
## 3 -0.09913143 -0.3463167 -0.05914688 -0.2027289  0.1457924       1
## 4 -0.48741352 -0.4638960 -0.47368490 -0.4738404 -0.3489628       1
## 5 -0.57744995 -0.4638960 -0.52603857 -0.4924883 -0.4196421       1
## 6 -0.34110433 -0.3306394 -0.37069409 -0.4135934 -0.4196421       1
# Visualización de los conglomerados
ggplot(cluster_complete, aes(x = HURTO, y = ROBO, color = cluster)) +
  geom_point(size = 2) +
  labs(title = "Conglomerados jerárquicos (3 clusters)", x = "Hurto (escalada)", y = "Robo (escalado)") +
  theme_minimal()

# Visualización con 3 clusters
fviz_cluster(list(data = variables_escaladas, cluster = cutree(hc, k = 3)), geom = "point", ellipse.type = "convex") +
  labs(title = "Distribución con 3 Clusters")

# Agrupación jerárquica con enlace simple
hc_single <- hclust(d, method = "single")

# Visualizamos el dendrograma
fviz_dend(hc_single, k = 3, rect = TRUE, show_labels = FALSE, cex = 0.6)

# Visualización de clusters con enlace simple
fviz_cluster(list(data = variables_escaladas, cluster = cutree(hc_single, k = 3)), geom = "point") +
  labs(title = "Clusters - Método Single")

#3.3. DBSCAN

DBSCAN son las siglas de Density-Based Spatial Clustering of Applications with Noise (agrupación espacial de aplicaciones con ruido basada en la densidad). Agrupa los puntos que están muy juntos y trata los puntos en regiones de baja densidad como valores atípicos.

library(dbscan)
## 
## Attaching package: 'dbscan'
## The following object is masked from 'package:stats':
## 
##     as.dendrogram

Realizamos el gráfico llamado “distancias kNN”. El eje x del gráfico representa los puntos de datos ordenados por su distancia al k-ésimo vecino más cercano (kNN), y el eje y muestra la distancia al k-ésimo vecino más cercano. En el contexto de DBSCAN, este tipo de gráfico puede ser útil para identificar un valor de epsilon (ε), que es un parámetro clave en DBSCAN que define la distancia máxima entre dos puntos para que uno sea considerado vecino del otro. Con el “abline” agregamos una línea horizontal al gráfico en la posición y = 0.32. Este valor está relacionado con la distancia epsilon (ε) que estamos considerando como umbral para definir vecinos. La línea roja podría ayudarte a visualizar y seleccionar un valor de epsilon que sea adecuado para tu análisis.

# Crear el gráfico de distancias kNN
kNNdistplot(variables_escaladas, k = 3)  # k = minPts - 1
abline(h = 0.32, col = "red")  # Valor sugerido para epsilon

Usamos los siguientes parámetros:

Valor de ε (epsilon): 0.32 Número mínimo de puntos (minPts): 4

minPts define cuántos puntos de la vecindad épsilon son necesarios para que un punto sea un punto central. A menudo se elige como parámetro el suavizado. Aquí utilizamos minPts = 4.

Para decidir el épsilon, se utiliza el gráfico de distancia kNN. Si la distancia entre dos puntos es menor o igual a 0.32, entonces estos puntos se consideran vecinos cercanos y se agruparán juntos en el mismo clúster si ambos tienen suficientes vecinos (determinado por el parámetro minPts). Si la distancia entre dos puntos es mayor que 0.32, entonces estos puntos no se consideran vecinos cercanos y pueden estar en diferentes clústeres, o si no tienen suficientes vecinos dentro de este radio, se considerarán puntos de ruido (outliers)

# Ejecutar DBSCAN
db <- dbscan(variables_escaladas, eps = 0.32, minPts = 4)

# Ver los resultados
db
## DBSCAN clustering for 271 objects.
## Parameters: eps = 0.32, minPts = 4
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 1 cluster(s) and 77 noise points.
## 
##   0   1 
##  77 194 
## 
## Available fields: cluster, eps, minPts, metric, borderPoints

Vemos la estructura de lo obtenido

# Agregamos los clusters a las variables escaladas
clustered_data <- as.data.frame(variables_escaladas) %>%
  add_column(cluster = factor(db$cluster))

# Visualización básica de los clusters
ggplot(clustered_data, aes(x = HURTO, y = ROBO, color = cluster)) +
  geom_point() +
  labs(title = "Resultados DBSCAN", x = "Hurto (escalada)", y = "Robo (escalado)") +
  theme_minimal()

library(factoextra)

# Visualizar los clusters con factoextra
fviz_cluster(db, variables_escaladas, geom = "point") +
  labs(title = "Visualización de Clusters DBSCAN")

3.4. Partición alrededor de los medoides (Partitioning Around Medoids) - PAM

# Seleccionamos una muestra aleatoria de 10,000 registros (ajustar según la memoria disponible)
set.seed(123)
#datos_muestra <- datos[sample(1:nrow(datos), 10000), ]

# Preparamos las variables mixtas (numéricas y categóricas)
data_limpia$Departamento <- as.factor(data_limpia$Departamento)

# Calculamos la matriz de distancias Gower
library(cluster)
d <- daisy(data_limpia[, c("ESTAFA", "EXTORSION", "HURTO", "ROBO", "HOMICIDIO", "Departamento")], metric = "gower")
# Aplicamos PAM
p <- pam(d, k = 3)

# Revisamos los resultados
p
## Medoids:
##       ID    
## [1,] 197 197
## [2,] 116 116
## [3,] 176 176
## Clustering vector:
##   [1] 1 1 2 1 1 1 1 1 3 1 1 1 2 1 1 1 3 2 2 3 1 1 2 2 2 2 1 2 1 2 1 1 1 1 1 2 1
##  [38] 1 1 1 1 2 1 1 1 1 1 2 1 1 3 2 2 3 1 2 1 2 2 2 1 1 1 1 1 1 1 1 2 1 1 2 1 1
##  [75] 1 1 3 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 2 1 2 1 1 1 1 1 1 2 1 1
## [149] 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## [186] 3 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 2 2 1 2 1 1 1 1
## [223] 3 3 1 1 1 1 2 1 2 1 1 1 1 2 3 1 2 1 3 1 3 1 2 1 1 1 1 1 2 1 1 1 1 2 1 1 2
## [260] 2 1 1 1 1 1 1 3 2 2 1 1
## Objective function:
##     build      swap 
## 0.1605865 0.1520545 
## 
## Available components:
## [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
## [6] "clusinfo"   "silinfo"    "diss"       "call"
# Añadimos los clusters al subconjunto
data_limpia$cluster <- factor(p$clustering)

# Vemos las primeras filas del subconjunto con los clusters
head(data_limpia)
## # A tibble: 6 × 10
##   `Ubigeo Hecho` Departamento Provincia   Distrito  ESTAFA EXTORSION HURTO  ROBO
##            <dbl> <fct>        <chr>       <chr>      <dbl>     <dbl> <dbl> <dbl>
## 1          10101 AMAZONAS     CHACHAPOYAS CHACHAPO…    127        34   172    31
## 2          10203 AMAZONAS     BAGUA       COPALLIN       2         1     8     5
## 3          10701 AMAZONAS     UTCUBAMBA   BAGUA GR…     86        16   552   203
## 4          10702 AMAZONAS     UTCUBAMBA   CAJARURO      17         1    69    14
## 5          20507 ÁNCASH       BOLOGNESI   COLQUIOC       1         1     8     1
## 6          20801 ÁNCASH       CASMA       CASMA         43        18   189    56
## # ℹ 2 more variables: HOMICIDIO <dbl>, cluster <fct>
library(ggplot2)

# Visualizamos en un gráfico de dispersión
ggplot(data_limpia, aes(x = HURTO, y = ROBO, color = cluster)) +
  geom_point(size = 2) +
  labs(title = "Clusters PAM con Variables Mixtas", x = "Hurto", y = "Robo") +
  theme_minimal()