Análisis de Clusters en Comercio entre EE.UU. y México usando DBSCAN

📂 1. Cargar datos

# Datos simulados para Alabama
df <- read_csv(here("data", "State Imports by NAICS Commodities.csv"), skip = 2) %>% 
  clean_names() %>% 
  select(-x6)

New names:
• `` -> `...6`

Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
  dat <- vroom(...)
  problems(dat)

Rows: 7992 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): State, Commodity, Country
dbl (1): Time
num (1): Customs Value (Gen) ($US)
lgl (1): ...6

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

🧹 2. Limpieza de datos

Responde a lo siguiente:

1. En este chunk borré los valores vacíos en lugar de convertirlos en 0, ¿por qué lo decidí de esta forma? HINT: Tiene que ver con el objetivo de este ejercicio “descubrir las agrupaciones naturales de commodities”.

2. En el siguiente chunk hago un muy elegante filtro usando str_detect, ¿por qué lo hice?

df_clean <- df %>%
  #drop_na(customs_value_gen_us) %>%
  mutate(customs_value_gen_us = ifelse(is.na(customs_value_gen_us), 0, customs_value_gen_us)) %>%
  filter(str_detect(commodity, "^\\d{3} ")) %>%  # Filtra exactamente 3 dígitos
  mutate(log_value = log1p(customs_value_gen_us)) %>% # Transformación logarítmica segura
  arrange(-customs_value_gen_us) %>% 
  head(1000)

📈 3. Exploración rápida

df_clean %>% 
  mutate(x_label = str_c(state, " - ", commodity)) %>% 
  slice_max(order_by = customs_value_gen_us, n = 10) %>% 
  ggplot(aes(x = fct_reorder(x_label, customs_value_gen_us), y = customs_value_gen_us)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  scale_y_continuous(labels = scales::comma) +
  labs(
    title = "Top 10 importaciones desde México por estado y producto (2024)",
    x = "Estado - Producto",
    y = "Valor en dólares"
  )

🧠 4. Aplicación de DBSCAN

DBSCAN agrupa por densidad sin necesidad de especificar número de clusters.

# Suponiendo que df_clean contiene tus datos de comercio
# Seleccionar múltiples variables relevantes para el agrupamiento (clustering)
data_matrix <- df_clean %>% 
  select(log_value, state, commodity) %>%
  # Convertir variables categóricas a numéricas si es necesario
  mutate(
    state_num = as.numeric(factor(state)),
    commodity_num = as.numeric(factor(commodity))
  ) %>%
  # Seleccionar las columnas numéricas para el agrupamiento
  select(log_value, state_num, commodity_num) %>%
  # Escalar los datos (importante para DBSCAN)
  scale()

# Probar diferentes valores de eps
# Primero, busquemos un eps apropiado usando un gráfico de distancia k-vecinos

kNNdistplot(data_matrix, k = 3)
abline(h = 0.6, col = "red") # Ajusta este valor según el "codo" del gráfico

# First create the proper data matrix for DBSCAN
# This step is crucial and was missing
data_matrix <- df_clean %>% 
  # Create numeric representations of categorical variables
  mutate(
    state_num = as.numeric(factor(state)),
    commodity_num = as.numeric(factor(commodity))
  ) %>%
  # Select features for clustering
  select(log_value, state_num, commodity_num) %>%
  # Scale the data for better clustering
  scale()

# Run DBSCAN with appropriate parameters
modelo_dbscan <- dbscan(data_matrix, eps = 0.5, minPts = 5)

# Add clusters to original data
df_clean <- df_clean %>%
  mutate(cluster = as.factor(modelo_dbscan$cluster))

# Analyze cluster results
table(df_clean$cluster)


  0   1   2   3 
 44 935  16   5

# For visualization, you need to use numeric values
# This plot will work correctly
ggplot(df_clean, aes(x = log_value, y = state, color = cluster)) +
  geom_point(alpha = 0.7) +
  theme_minimal() +
  labs(title = "DBSCAN Clustering of Mexico-US Trade Data")

# If you want to see commodity patterns clearly, try this visualization
ggplot(df_clean, aes(x = log_value, y = commodity, color = cluster)) +
  geom_point(alpha = 0.7) +
  theme_minimal() +
  labs(title = "Commodity Clusters in Mexico-US Trade Data")

# This analysis to examine cluster characteristics is good
cluster_summary <- df_clean %>%
  group_by(cluster) %>%
  summarize(
    count = n(),
    avg_value = mean(exp(log_value)),
    top_commodities = paste(names(sort(table(commodity), decreasing = TRUE)[1:5]), collapse = ", "),
    top_states = paste(names(sort(table(state), decreasing = TRUE)[1:5]), collapse = ", ")
  )

knitr::kable(cluster_summary)

cluster	count	avg_value	top_commodities	top_states
0	44	6030169619	111 Agricultural Products, 311 Food & Kindred Products, 336 Transportation Equipment, 112 Livestock & Livestock Products, 334 Computer & Electronic Products	Texas, California, Michigan, Georgia, Illinois
1	935	243306985	333 Machinery, Except Electrical, 326 Plastics & Rubber Products, 334 Computer & Electronic Products, 335 Electrical Equipment, Appliances & Components, 332 Fabricated Metal Products, Nesoi	Pennsylvania, New York, Ohio, Arizona, Illinois
2	16	573930418	312 Beverages & Tobacco Products, 311 Food & Kindred Products, 211 Oil & Gas, 114 Fish, Fresh/chilled/frozen & Other Marine Products, 112 Livestock & Livestock Products	California, Florida, Arizona, Georgia, Alabama
3	5	744497772	311 Food & Kindred Products, 312 Beverages & Tobacco Products, 211 Oil & Gas, NA, NA	New Jersey, New York, Pennsylvania, NA, NA

# Add this to see how commodities are distributed across clusters
commodity_distribution <- df_clean %>%
  group_by(commodity, cluster) %>%
  summarize(
    count = n(),
    avg_value = mean(exp(log_value))
  ) %>%
  arrange(commodity, cluster)

`summarise()` has grouped output by 'commodity'. You can override using the
`.groups` argument.

print(commodity_distribution)

# A tibble: 54 × 4
# Groups:   commodity [31]
   commodity                                             cluster count avg_value
   <chr>                                                 <fct>   <int>     <dbl>
 1 111 Agricultural Products                             0          16    1.19e9
 2 111 Agricultural Products                             1          18    5.85e7
 3 112 Livestock & Livestock Products                    0           3    2.79e8
 4 112 Livestock & Livestock Products                    1           4    3.82e7
 5 112 Livestock & Livestock Products                    2           1    4.27e8
 6 113 Forestry Products, Nesoi                          0           1    2.10e7
 7 113 Forestry Products, Nesoi                          1           3    5.56e6
 8 114 Fish, Fresh/chilled/frozen & Other Marine Produc… 0           1    1.88e6
 9 114 Fish, Fresh/chilled/frozen & Other Marine Produc… 1          11    7.14e6
10 114 Fish, Fresh/chilled/frozen & Other Marine Produc… 2           2    1.56e8
# ℹ 44 more rows

🎯 5. Visualización de clusters

library(umap)

# Aplicar UMAP para reducir dimensionalidad
set.seed(123) # Para reproducibilidad
umap_result <- umap(data_matrix, n_neighbors = 15, min_dist = 0.1, n_components = 2)

# Crear un dataframe con las coordenadas de UMAP y la información de clusters
umap_df <- data.frame(
  UMAP1 = umap_result$layout[,1],
  UMAP2 = umap_result$layout[,2]
)

# Añadir información de clusters y categorías
umap_df$cluster <- df_clean$cluster
umap_df$commodity <- df_clean$commodity
umap_df$state <- df_clean$state
umap_df$log_value <- df_clean$log_value

# Visualizar los clusters usando UMAP
ggplot(umap_df, aes(x = UMAP1, y = UMAP2, color = cluster)) +
  geom_point(alpha = 0.7, size = 2) +
  theme_minimal() +
  labs(title = "UMAP Visualization of DBSCAN Clusters in Mexico-US Trade Data",
       color = "Cluster") +
  scale_color_brewer(palette = "Set1")

# Visualización alternativa: colorear por commodity
# Limitar a las 10 commodities más frecuentes para claridad
top_commodities <- names(sort(table(df_clean$commodity), decreasing = TRUE)[1:10])
umap_df_top <- umap_df %>% filter(commodity %in% top_commodities)

ggplot(umap_df_top, aes(x = UMAP1, y = UMAP2, color = commodity)) +
  geom_point(alpha = 0.7, size = 2) +
  theme_minimal() +
  labs(title = "UMAP Visualization by Commodity (Top 10)",
       color = "Commodity") +
  theme(legend.position = "right")

# También puedes visualizar una versión que combine clusters y log_value
ggplot(umap_df, aes(x = UMAP1, y = UMAP2, color = cluster, size = log_value)) +
  geom_point(alpha = 0.6) +
  theme_minimal() +
  labs(title = "UMAP Visualization by Cluster and Trade Value",
       color = "Cluster", 
       size = "Log Value") +
  scale_color_brewer(palette = "Set1")

📌 6. Conclusiones rápidas

DBSCAN detecta agrupaciones sin definir el número de grupos.
Algunos productos quedaron como ruido (cluster 0).

🔥 7. Ejercicio adicional: Encontrar Estados Similares

Supongamos que tenemos datos de varios estados (simularemos aquí para el ejemplo).
Queremos encontrar estados que importan productos similares en valor.

df_clean <- df %>%
  mutate(customs_value_gen_us = ifelse(is.na(customs_value_gen_us), 0, customs_value_gen_us)) %>%
  filter(str_detect(commodity, "^\\d{4} ")) %>%  # Filtra exactamente 4 dígitos
  mutate(log_value = scale(customs_value_gen_us)) %>% # Transformación logarítmica segura
  arrange(-customs_value_gen_us) %>% 
  head(1000)

# Crear el dataframe 'df_agg' como antes
df_agg <- df_clean %>%
  group_by(state, commodity) %>%
  summarise(total_value = sum(customs_value_gen_us, na.rm = TRUE), .groups = "drop")

# Crear el dataframe 'df_wide' para pivotar las columnas
df_wide <- df_agg %>%
  pivot_wider(names_from = commodity, values_from = total_value, values_fill = 0)

# Preparamos la matriz para clustering
state_matrix <- df_wide %>% select(-state) %>% as.matrix()

# Asignar los nombres de los estados a las filas de la matriz
rownames(state_matrix) <- df_wide$state

# Cálculo de distancias
state_dist <- dist(state_matrix, method = "euclidean")

# Clustering jerárquico
hc <- hclust(state_dist, method = "ward.D2")

# Dendrograma con nombres de los estados como etiquetas
plot(hc, main = "Estados agrupados por patrón de importaciones", xlab = "", sub = "", labels = rownames(state_matrix))

# Aplicar transformación logarítmica para manejar valores de comercio sesgados
state_matrix <- df_wide %>% 
  select(-state) %>%
  select(where(~ sd(.) > 0)) %>%
  as.matrix()

# Agregar una constante pequeña antes del logaritmo para manejar ceros
state_matrix_log <- log1p(state_matrix)  # log(1+x) maneja ceros de forma segura
rownames(state_matrix_log) <- df_wide$state

# 2. Escalar los datos para el agrupamiento (clustering)
state_matrix_scaled <- scale(state_matrix_log)

# 3. Probar diferentes valores de eps para encontrar el agrupamiento óptimo
# Optimizar el parámetro eps usando el gráfico de distancia k-vecinos
kNNdistplot(state_matrix_scaled, k = 3)
abline(h = 15, col = "red")  # Elegir eps en el punto de “codo” del gráfico

# 4. Ejecutar DBSCAN con parámetros optimizados
dbscan_result <- dbscan(state_matrix_scaled, eps = 15, minPts = 5)

# 5. Crear múltiples visualizaciones para obtener mejores perspectivas

# A. Visualización con PCA (versión mejorada de la original)
df_pca <- data.frame(
  prcomp(state_matrix_scaled)$x[, 1:2],
  Cluster = as.factor(dbscan_result$cluster),
  State = df_wide$state
)

# Contar los estados por clúster para información del título
cluster_counts <- table(dbscan_result$cluster)
title_info <- paste("Agrupación de Estados por Patrón Comercial (DBSCAN)\nClústeres:", 
                   length(cluster_counts)-1, 
                   "| Puntos de ruido:", 
                   cluster_counts["0"])

# Gráfico PCA mejorado con estética optimizada
pca_plot <- ggplot(df_pca, aes(x = PC1, y = PC2, color = Cluster, label = State)) +
  geom_point(size = 3, alpha = 0.7) +
  ggrepel::geom_text_repel(
    size = 3.5, 
    box.padding = 0.35, 
    point.padding = 0.5,
    segment.color = "grey50"
  ) +
  theme_minimal() +
  labs(
    title = title_info,
    subtitle = "Análisis de Componentes Principales (PCA)",
    x = "Componente Principal 1",
    y = "Componente Principal 2"
  ) +
  scale_color_viridis_d(option = "plasma", end = 0.9) +
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(face = "italic")
  )

print(pca_plot)

Warning: ggrepel: 39 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

# B. Visualización con UMAP (a menudo mejor para preservar estructuras locales)
set.seed(42)  # Para reproducibilidad
umap_result <- umap(state_matrix_scaled, n_neighbors = 5, min_dist = 0.1, n_components = 2)

df_umap <- data.frame(
  UMAP1 = umap_result$layout[,1],
  UMAP2 = umap_result$layout[,2],
  Cluster = as.factor(dbscan_result$cluster),
  State = df_wide$state
)

umap_plot <- ggplot(df_umap, aes(x = UMAP1, y = UMAP2, color = Cluster, label = State)) +
  geom_point(size = 3, alpha = 0.7) +
  ggrepel::geom_text_repel(
    size = 3.5, 
    box.padding = 0.35, 
    point.padding = 0.5,
    segment.color = "grey50"
  ) +
  theme_minimal() +
  labs(
    title = title_info,
    subtitle = "Reducción de Dimensionalidad con UMAP",
    x = "Dimensión UMAP 1",
    y = "Dimensión UMAP 2"
  ) +
  scale_color_viridis_d(option = "plasma", end = 0.9) +
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold"),
    plot.subtitle = element_text(face = "italic")
  ) 

print(umap_plot)

🧠 Interpretación

Estados cercanos en el dendrograma tienen patrones de importaciones similares.
Esto puede ayudar a diseñar políticas regionales o detectar estrategias comerciales compartidas.

✅ Resumen final

DBSCAN encuentra agrupamientos de productos sin definir el número de clusters.
Podemos usar clustering jerárquico para encontrar estados con patrones similares de importaciones.
El análisis de densidades y patrones de comercio revela oportunidades y riesgos en relaciones internacionales.