USArrests

Contexto

La base de datos USArrests contiene estadísticas en arrestos por cada 100,000 residentes por agresión, asesinato y violación en cada uno de los 50 estados de EE.UU. en 1973.

Instalar paquetes y llamar librerías

#install.packages("cluster") #Para agrupamiento de datos 
library(cluster)
#install.packages("ggplot2") #Gráficas
library(ggplot2)
#install.packages("factoextra") #Visualizar clusters
library(factoextra)
#install.packages("data.table") #Conjunto de datos grande 
library(data.table)
#install.packages("tidyverse")
library(tidyverse)
#install.packages("caret") # Algoritmos de aprendizaje automático
library(caret)
#install.packages("datasets") # Para usar base de datos "iris"
library(datasets)
#install.packages("ggplot2") # Graficas con mejor diseño
library(ggplot2)
#install.packages("lattice") # Crear gráficos
library(lattice)
#install.packages("DataExplorer") # Análisis Descriptivo
library(DataExplorer)
#install.packages("kernlab")
library(kernlab)
#install.packages("randomForest")
library(randomForest)
library(tibble)
#install.packages("sf") #Analisis de datos espaciales
library(sf)
#install.packages("rnaturalearth") #Proporciona límites geograficos
library(rnaturalearth)
#install.packages("rnaturalearthdata") #Datos de geografia
library(rnaturalearthdata)
#install.packages("devtools")
library(devtools)
devtools::install_github("ropensci/rnaturalearthhires") #Mapa de México particular

Importar base de datos

df <- USArrests

Análisis descriptivo

create_report(df)

##   |                                             |                                     |   0%  |                                             |.                                    |   2%                                   |                                             |..                                   |   5% [global_options]                  |                                             |...                                  |   7%                                   |                                             |....                                 |  10% [introduce]                       |                                             |....                                 |  12%                                   |                                             |.....                                |  14% [plot_intro]

##   |                                             |......                               |  17%                                   |                                             |.......                              |  19% [data_structure]                  |                                             |........                             |  21%                                   |                                             |.........                            |  24% [missing_profile]

##   |                                             |..........                           |  26%                                   |                                             |...........                          |  29% [univariate_distribution_header]  |                                             |...........                          |  31%                                   |                                             |............                         |  33% [plot_histogram]

##   |                                             |.............                        |  36%                                   |                                             |..............                       |  38% [plot_density]                    |                                             |...............                      |  40%                                   |                                             |................                     |  43% [plot_frequency_bar]              |                                             |.................                    |  45%                                   |                                             |..................                   |  48% [plot_response_bar]               |                                             |..................                   |  50%                                   |                                             |...................                  |  52% [plot_with_bar]                   |                                             |....................                 |  55%                                   |                                             |.....................                |  57% [plot_normal_qq]

##   |                                             |......................               |  60%                                   |                                             |.......................              |  62% [plot_response_qq]                |                                             |........................             |  64%                                   |                                             |.........................            |  67% [plot_by_qq]                      |                                             |..........................           |  69%                                   |                                             |..........................           |  71% [correlation_analysis]

##   |                                             |...........................          |  74%                                   |                                             |............................         |  76% [principal_component_analysis]

##   |                                             |.............................        |  79%                                   |                                             |..............................       |  81% [bivariate_distribution_header]   |                                             |...............................      |  83%                                   |                                             |................................     |  86% [plot_response_boxplot]           |                                             |.................................    |  88%                                   |                                             |.................................    |  90% [plot_by_boxplot]                 |                                             |..................................   |  93%                                   |                                             |...................................  |  95% [plot_response_scatterplot]       |                                             |.................................... |  98%                                   |                                             |.....................................| 100% [plot_by_scatterplot]           

## /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/pandoc +RTS -K512m -RTS /Users/mariadelbosque/Desktop/TEC/CONCENTRACION/R/report.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /Users/mariadelbosque/Desktop/TEC/CONCENTRACION/R/report.html --lua-filter /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/rmarkdown/rmarkdown/lua/latex-div.lua --embed-resources --standalone --variable bs3=TRUE --section-divs --table-of-contents --toc-depth 6 --template /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=yeti --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /var/folders/fg/qsq_j5fx5wz0c1xt7brrw6fh0000gn/T//Rtmpa85o7o/rmarkdown-str36e61de02cfb.html

plot_missing(df)

plot_histogram(df)

plot_correlation(df)

# Clusters

Escalar la base de datos

df_escalado <-scale(df)
df_escalado <- subset(df_escalado,select = -UrbanPop)
df_escalado <-scale(df_escalado)

Generar los segmentos

grupos <- 4 #Inicio con cualquier valor y luego verifico
segmentos <- kmeans(df_escalado, grupos)

Asignar grupos a los datos

asignacion <- cbind(df, cluster=segmentos$cluster)

Graficar los clusters

fviz_cluster(segmentos, data = df)

## Optimizar la cantidad de grupos

#La cantidad óptima de grupos corresponde al primer punto más alto de la gráfica
set.seed(123)
#optimizacion <- clusGap(df_escalado, FUN=kmeans, nstart=1, K.max =10)
#plot(optimizacion, xlab = "Número de Clusters k")

# Crear una función para calcular la suma de los cuadrados dentro del grupo (WCSS)
wcss <- function(k) {
  kmeans(df_escalado, centers = k, nstart = 10)$tot.withinss
}

# Evaluar WCSS para diferentes valores de k
k_values <- 1:10
wcss_values <- sapply(k_values, wcss)

# Graficar el método del codo
plot(k_values, wcss_values, type = "b", pch = 19, frame = FALSE,
     xlab = "Número de Clusters k",
     ylab = "WCSS (Suma de los Cuadrados Dentro del Grupo)",
     main = "Método del Codo para Optimización de Clusters")

## Comparar segmentos

asignacion$state <- rownames(asignacion)

# 4. Calculate average metrics by cluster for ranking
cluster_avg <- asignacion %>%
  group_by(cluster) %>%
  summarise(
    Murder = mean(Murder, na.rm = TRUE),
    Assault = mean(Assault, na.rm = TRUE),
    Rape = mean(Rape, na.rm = TRUE)
  ) %>%
  mutate(
    Suma = Murder + Assault + Rape,
    Ranking = rank(-Suma, ties.method = "min")  # Higher crime gets lower ranking (1)
  )

# 5. Create security level labels based on ranking
cluster_labels <- cluster_avg %>%
  mutate(nivel_seguridad = case_when(
    Ranking == 1 ~ "Muy inseguro",
    Ranking == 2 ~ "Inseguro", 
    Ranking == 3 ~ "Seguro",
    Ranking == 4 ~ "Muy seguro",
    TRUE ~ "Desconocido"
  )) %>%
  select(cluster, Ranking, nivel_seguridad)

# 6. Join labels back to the data
asignacion <- left_join(asignacion, cluster_labels, by = "cluster")

# 7. Get US map
us_map <- ne_states(country = "United States of America", returnclass = "sf")

# 8. Create manual mapping between map state names and USArrests state names
state_mapping <- data.frame(
  map_name = c(
    "Alabama", "Alaska", "Arizona", "Arkansas", "California", 
    "Colorado", "Connecticut", "Delaware", "Florida", "Georgia", 
    "Hawaii", "Idaho", "Illinois", "Indiana", "Iowa", 
    "Kansas", "Kentucky", "Louisiana", "Maine", "Maryland", 
    "Massachusetts", "Michigan", "Minnesota", "Mississippi", "Missouri", 
    "Montana", "Nebraska", "Nevada", "New Hampshire", "New Jersey", 
    "New Mexico", "New York", "North Carolina", "North Dakota", "Ohio", 
    "Oklahoma", "Oregon", "Pennsylvania", "Rhode Island", "South Carolina", 
    "South Dakota", "Tennessee", "Texas", "Utah", "Vermont", 
    "Virginia", "Washington", "West Virginia", "Wisconsin", "Wyoming"
  ),
  data_name = state.name  # This matches the row names in USArrests
)

# 9. Add mapping column to the map data
us_map$data_name <- NA
for (i in 1:nrow(state_mapping)) {
  idx <- which(us_map$name == state_mapping$map_name[i])
  if (length(idx) > 0) {
    us_map$data_name[idx] <- state_mapping$data_name[i]
  }
}

# 10. Join map with data
us_clustered <- left_join(us_map, asignacion, by = c("data_name" = "state"))

# 11. Define colors for the security levels
nivel_colors <- c(
  "Muy inseguro" = "red", 
  "Inseguro" = "orange", 
  "Seguro" = "yellow", 
  "Muy seguro" = "darkgreen",
  "Desconocido" = "gray"
)

# 12. Print summary to check
print("Cluster summaries with rankings:")

## [1] "Cluster summaries with rankings:"

print(cluster_labels)

## # A tibble: 4 × 3
##   cluster Ranking nivel_seguridad
##     <int>   <int> <chr>          
## 1       1       1 Muy inseguro   
## 2       2       2 Inseguro       
## 3       3       4 Muy seguro     
## 4       4       3 Seguro

print("Distribution of security levels:")

## [1] "Distribution of security levels:"

print(table(asignacion$nivel_seguridad))

## 
##     Inseguro Muy inseguro   Muy seguro       Seguro 
##           12           19            7           12

# 13. Create the map
ggplot(data = us_clustered) +
  geom_sf(aes(fill = nivel_seguridad), color = "black", size = 0.2) +
  scale_fill_manual(values = nivel_colors, 
                    name = "Nivel de Seguridad",
                    na.value = "gray") +
  labs(title = "Mapa de Seguridad en EE.UU. (1973)",
       subtitle = "Clasificación basada en tasas de criminalidad",
       caption = "Fuente: USArrests") +
  theme_minimal()

** NOTA: La variable que queremos predecir debe tener formato de FACTOR - en este caso serían los clusters**

CARET

Partir los datos 80-20

set.seed(123)
renglones_entrenamiento <- createDataPartition(asignacion$cluster, p=0.8, list = FALSE)
entrenamiento <- asignacion[renglones_entrenamiento, ] 
prueba <- asignacion[-renglones_entrenamiento, ]

entrenamiento <- entrenamiento %>% select(-state)
prueba <- prueba %>% select(-state)

Modelo 1. SVM Linear

# Asegurar que la variable dependiente 'cluster' es un factor en entrenamiento y prueba
entrenamiento$cluster <- factor(entrenamiento$cluster)
prueba$cluster <- factor(prueba$cluster)


# Entrenar el modelo SVM
modelo1 <- train(cluster ~ ., data = entrenamiento,  
                 method = "svmLinear", 
                 preProcess = c("scale", "center"),
                 trControl = trainControl(method = "cv", number = 10),
                 tuneGrid = data.frame(C = 1)  
)

# Hacer predicciones
resultado_entrenamiento1 <- predict(modelo1, entrenamiento)
resultado_prueba1 <- predict(modelo1, prueba)

# Convertir las predicciones en factor con los mismos niveles que 'entrenamiento$cluster'
resultado_entrenamiento1 <- factor(resultado_entrenamiento1, levels = levels(entrenamiento$cluster))
resultado_prueba1 <- factor(resultado_prueba1, levels = levels(prueba$cluster))

# Matriz de Confusión del Entrenamiento
mcre1 <- confusionMatrix(resultado_entrenamiento1, entrenamiento$cluster)
print(mcre1)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  1  2  3  4
##          1 15  0  0  0
##          2  0 10  0  0
##          3  0  0  6  0
##          4  0  0  0 10
## 
## Overall Statistics
##                                     
##                Accuracy : 1         
##                  95% CI : (0.914, 1)
##     No Information Rate : 0.3659    
##     P-Value [Acc > NIR] : < 2.2e-16 
##                                     
##                   Kappa : 1         
##                                     
##  Mcnemar's Test P-Value : NA        
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3 Class: 4
## Sensitivity            1.0000   1.0000   1.0000   1.0000
## Specificity            1.0000   1.0000   1.0000   1.0000
## Pos Pred Value         1.0000   1.0000   1.0000   1.0000
## Neg Pred Value         1.0000   1.0000   1.0000   1.0000
## Prevalence             0.3659   0.2439   0.1463   0.2439
## Detection Rate         0.3659   0.2439   0.1463   0.2439
## Detection Prevalence   0.3659   0.2439   0.1463   0.2439
## Balanced Accuracy      1.0000   1.0000   1.0000   1.0000

# Matriz de Confusión del Resultado de la Prueba
mcrp1 <- confusionMatrix(resultado_prueba1, prueba$cluster)
print(mcrp1)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 1 2 3 4
##          1 4 0 0 0
##          2 0 2 0 0
##          3 0 0 1 0
##          4 0 0 0 2
## 
## Overall Statistics
##                                      
##                Accuracy : 1          
##                  95% CI : (0.6637, 1)
##     No Information Rate : 0.4444     
##     P-Value [Acc > NIR] : 0.0006766  
##                                      
##                   Kappa : 1          
##                                      
##  Mcnemar's Test P-Value : NA         
## 
## Statistics by Class:
## 
##                      Class: 1 Class: 2 Class: 3 Class: 4
## Sensitivity            1.0000   1.0000   1.0000   1.0000
## Specificity            1.0000   1.0000   1.0000   1.0000
## Pos Pred Value         1.0000   1.0000   1.0000   1.0000
## Neg Pred Value         1.0000   1.0000   1.0000   1.0000
## Prevalence             0.4444   0.2222   0.1111   0.2222
## Detection Rate         0.4444   0.2222   0.1111   0.2222
## Detection Prevalence   0.4444   0.2222   0.1111   0.2222
## Balanced Accuracy      1.0000   1.0000   1.0000   1.0000

Resumen de Resultados

resultados <- data.frame(
  "SVM Lineal" = c(mcre1$overall["Accuracy"], mcrp1$overall["Accuracy"])
)
rownames(resultados) <- c("Precision de Entrenamiento", "Precision de Prueba")
resultados

##                            SVM.Lineal
## Precision de Entrenamiento          1
## Precision de Prueba                 1