Examen:

Datos. Utilice la librería wbstats para descargar, desde la API del Banco Mundial, los siguientes indicadores para el periodo 2015-2024:

Crecimiento del PIB (% anual): NY.GDP.MKTP.KD.ZG Desempleo (% de la fuerza laboral): SL.UEM.TOTL.ZS Remesas personales recibidas (% del PIB): BX.TRF.PWKR.DT.GD.ZS Exportaciones de bienes y servicios (% del PIB): NE.EXP.GNFS.ZS

Pregunta 1. Preparación de los datos

Seleccione un grupo de 20 países de su elección y descargue los indicadores indicados.
Construya la matriz de datos tomando, para cada país e indicador, el último dato disponible.
Trate los valores faltantes sustituyéndolos por la media de la columna, y verifique que no queden valores NA.
Estandarice las variables mediante scale().

Pregunta 2. Número óptimo de clústeres

Sobre la matriz estandarizada obtenida en la pregunta 1, determine el número óptimo de clústeres (k) para cada método, usando el criterio indicado:

K-Means: método del codo.
K-Medoides (PAM): coeficiente de silueta.
Ward: criterio del salto de altura en el dendrograma.

Presente el gráfico de cada criterio.

Pregunta 3. Aplicación de los métodos con el k seleccionado

Sobre la matriz estandarizada de la pregunta 1, aplique los siguientes métodos de agrupamiento (DBSCAN se aplica sobre las variables de comercio exterior indicadas en el inciso d):

K-Means, usando el k obtenido con el método del codo en la pregunta 2.
K-Medoides (PAM), usando el k obtenido con el coeficiente de silueta en la pregunta 2.
Ward, usando el k obtenido con el criterio del salto de altura en la pregunta 2.
DBSCAN, usando las variables Exportaciones_USD (NE.EXP.GNFS.CD) e Importaciones_USD (NE.IMP.GNFS.CD). Como este método no requiere un k previo, determine sus parámetros minPts y epsilon.

Pregunta 4. Representación gráfica de los clústeres

Grafique los clústeres obtenidos en la pregunta 3 con los cuatro métodos (K-Means, K-Medoides PAM, DBSCAN y Ward).

Datos

library(wbstats)
library(dplyr)
library(tidyr)
library(tibble)
library(knitr)
library(kableExtra)
library(factoextra)
library(cluster)
library(dbscan)
library(ggplot2)
library(ggrepel)

Parte 1

a) Selección de 20 países

paises_wb <- c(
  # América Latina
  "ARG","BRA","CHL","COL","CRI","ECU","MEX","PER","SLV","URY",
  # Norteamérica y Asia-Pacífico
  "USA","CAN","JPN","KOR",
  # Europa
  "ESP","FRA","DEU","ITA","GBR","POL"
)

Descarga de indicadores

indicadores <- c(
  Crecimiento_PIB = "NY.GDP.MKTP.KD.ZG",
  Desempleo       = "SL.UEM.TOTL.ZS",
  Remesas_personales     = "BX.TRF.PWKR.DT.GD.ZS",
  Exportaciones   = "NE.EXP.GNFS.ZS"
)

datos_wb <- wb_data(
  indicator   = indicadores,
  country     = paises_wb,
  start_date  = 2015,
  end_date    = 2024,
  return_wide = TRUE
)

b) matriz

ultimo_dato <- datos_wb %>%
  select(iso3c, country, date, all_of(names(indicadores))) %>%
  pivot_longer(cols = all_of(names(indicadores)),
               names_to = "Indicador", values_to = "Valor") %>%
  filter(!is.na(Valor)) %>%
  group_by(iso3c, country, Indicador) %>%
  slice_max(date, n = 1) %>%
  ungroup()

data_cluster <- ultimo_dato %>%
  select(iso3c, country, Indicador, Valor) %>%
  pivot_wider(names_from = Indicador, values_from = Valor)

data_cluster %>%
  kable(caption = "Matriz de variables por país (último dato disponible)") %>%
  kable_classic(html_font = "Times New Roman", font_size = 13)

Matriz de variables por país (último dato disponible)
iso3c	country	Crecimiento_PIB	Desempleo	Exportaciones	Remesas_personales
ARG	Argentina	-1.3429307	7.150	15.24975	0.1635999
BRA	Brazil	3.4193152	6.801	17.93858	0.2242818
CAN	Canada	2.0462661	6.351	32.46912	0.0375070
CHL	Chile	2.8051264	8.718	33.73207	0.0304727
COL	Colombia	1.4933073	9.619	16.08045	2.8235717
CRI	Costa Rica	4.0829244	6.936	37.52984	0.7494093
DEU	Germany	-0.4958519	3.400	41.43406	0.4731113
ECU	Ecuador	-1.9440902	3.453	30.49788	5.2861312
ESP	Spain	3.4552541	11.400	37.05588	0.3667485
FRA	France	1.1904680	7.400	33.88975	1.2268965
GBR	United Kingdom	1.0802772	4.361	31.05614	0.1308153
ITA	Italy	0.7830391	6.500	32.33507	0.5094347
JPN	Japan	-0.2400842	2.500	21.97991	0.1108693
KOR	Korea, Rep.	2.0036110	2.784	44.35824	0.3883463
MEX	Mexico	1.3505763	2.678	37.28362	3.6950729
PER	Peru	3.5152749	5.199	28.88393	1.6910867
POL	Poland	3.0284163	2.807	52.17981	0.9468633
SLV	El Salvador	2.5904637	3.261	32.65986	24.3362061
URY	Uruguay	3.3257591	8.209	28.23490	0.1653514
USA	United States	2.7931872	4.022	10.97470	0.0297358

c) Tratamiento de datos faltantes (imputación por la media)

vars <- names(indicadores)

data_cluster_imputado <- data_cluster %>%
  mutate(across(all_of(vars),
                ~ ifelse(is.na(.x), mean(.x, na.rm = TRUE), .x)))

# Verificación: no deben quedar NA
data_cluster_imputado %>%
  summarise(across(all_of(vars), ~ sum(is.na(.x))))

## # A tibble: 1 × 4
##   Crecimiento_PIB Desempleo Remesas_personales Exportaciones
##             <int>     <int>              <int>         <int>
## 1               0         0                  0             0

d) Estandarización

matriz_paises <- data_cluster_imputado %>%
  select(-country) %>%
  column_to_rownames(var = "iso3c")

matriz_tipificada <- scale(matriz_paises)

Parte 2

a) K-Means: método del codo.

set.seed(123)

fviz_nbclust(matriz_tipificada, kmeans, method = "wss", k.max = 10) +
  labs(title = "Método del Codo para K-Means",
       subtitle = "Crecimiento del PIB, Desempleo, Remesas personales y Exportaciones",
       x = "Número de clústeres (k)", y = "WCSS (Suma de cuadrados intra-clúster)") +
  theme_minimal()

b) K-Medoides (PAM): coeficiente de silueta.

fviz_nbclust(matriz_tipificada, pam, method = "silhouette", k.max = 10,
             diss = dist(matriz_tipificada, method = "manhattan")) +
  labs(title = "Coeficiente de Silueta promedio por k (PAM, distancia Manhattan)") +
  theme_minimal()

c) Ward: criterio del salto de altura en el dendrograma.

distancias  <- dist(matriz_tipificada, method = "euclidean")
modelo_ward <- hclust(distancias, method = "ward.D2")

alturas      <- rev(modelo_ward$height)
saltos       <- -diff(alturas)
n_candidatos <- min(10, length(alturas))

plot(1:n_candidatos, alturas[1:n_candidatos], type = "b", pch = 19, col = "steelblue",
     xlab = "Número de clústeres", ylab = "Altura de fusión",
     main = "Identificación del salto más grande - Ward")

k_ward <- which.max(saltos[1:(n_candidatos - 1)]) + 1
abline(v = k_ward, col = "red", lty = 2)
text(k_ward, max(alturas[1:n_candidatos]) * 0.9,
     labels = paste0("k óptimo ≈ ", k_ward), col = "red", pos = 4)

k_ward

## [1] 3

Parte 3

a) Aplicación de K-Means con el k seleccionado

k_optimo <-5

set.seed(123)
resultado_kmeans <- kmeans(matriz_tipificada, centers = k_optimo, nstart = 25)

resultado_kmeans

## K-means clustering with 5 clusters of sizes 5, 3, 3, 8, 1
## 
## Cluster means:
##   Crecimiento_PIB  Desempleo Exportaciones Remesas_personales
## 1      -0.2067681 -0.9469142     1.0174696        -0.19297468
## 2       0.4804115  0.4354591    -1.5346069        -0.21166780
## 3      -1.7090163 -0.5018319    -0.7982709        -0.05844998
## 4       0.5283076  0.7324413     0.2162136        -0.29103768
## 5       0.4931941 -0.9258414     0.1815766         4.10352822
## 
## Clustering vector:
## ARG BRA CAN CHL COL CRI DEU ECU ESP FRA GBR ITA JPN KOR MEX PER POL SLV URY USA 
##   3   2   4   4   2   4   1   3   4   4   1   4   3   1   1   4   1   5   4   2 
## 
## Within cluster sum of squares by cluster:
## [1] 5.231273 3.371699 3.990399 7.898235 0.000000
##  (between_SS / total_SS =  73.0 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

b) Aplicación de K-Medoides (PAM) con el k seleccionado

k_op <- 2

set.seed(123)
modelo_pam <- pam(matriz_tipificada, k = k_op, metric = "manhattan")

modelo_pam

## Medoids:
##     ID Crecimiento_PIB  Desempleo Exportaciones Remesas_personales
## CAN  3       0.1749824  0.2580647     0.1630425         -0.3946317
## MEX 15      -0.2318120 -1.1492127     0.6308600          0.2824546
## Clustering vector:
## ARG BRA CAN CHL COL CRI DEU ECU ESP FRA GBR ITA JPN KOR MEX PER POL SLV URY USA 
##   1   1   1   1   1   1   2   2   1   1   1   1   2   2   2   1   2   2   1   1 
## Objective function:
##    build     swap 
## 2.319271 2.319271 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"

c) Aplicación de Ward con el k seleccionado

grupos_ward <- cutree(modelo_ward, k = k_ward)
table(grupos_ward)

## grupos_ward
##  1  2  3 
## 12  7  1

data.frame(
  k_clusteres     = 1:n_candidatos,
  altura_fusion   = round(alturas[1:n_candidatos], 2),
  salto_siguiente = c(round(saltos[1:(n_candidatos - 1)], 2), NA)
) %>%
  kable(
    col.names = c("Fusión j (pasa de j + 1 a j grupos)", "Altura de fusión (ESS)", "Salto a la fusión siguiente"),
    caption = "Alturas de fusión y saltos - Ward"
  ) %>%
  kable_classic(html_font = "Times New Roman", font_size = 14) %>%
  kable_styling()

Alturas de fusión y saltos - Ward
Fusión j (pasa de j + 1 a j grupos)	Altura de fusión (ESS)	Salto a la fusión siguiente
1	6.61	0.82
2	5.79	1.47
3	4.32	0.65
4	3.67	0.66
5	3.01	0.43
6	2.58	0.31
7	2.26	0.17
8	2.10	0.12
9	1.97	0.32
10	1.65	NA

d) Aplicación de DBSCAN

DBSCAN agrupa por densidad, por lo que funciona mejor con pocas dimensiones. Se utilizan dos variables derivadas del comercio exterior, construidas con las exportaciones (NE.EXP.GNFS.CD) y las importaciones (NE.IMP.GNFS.CD) de bienes y servicios en dólares corrientes:

log_total = ln(X + M + 1): tamaño del comercio exterior del país.
log_ratio = ln(X + 1) - ln(M + 1): balanza comercial relativa (positiva en exportadores netos, negativa en importadores netos).

Se aplican logaritmos porque el comercio en dólares tiene un fuerte sesgo a la derecha.

Descarga y preparación de los datos

ind_db <- c(
  Exportaciones_USD = "NE.EXP.GNFS.CD",
  Importaciones_USD = "NE.IMP.GNFS.CD"
)

datos_comercio <- wb_data(
  indicator   = ind_db,
  country     = paises_wb,
  start_date  = 2015,
  end_date    = 2024,
  return_wide = TRUE
)

datos_db <- datos_comercio %>%
  filter(!is.na(Exportaciones_USD), !is.na(Importaciones_USD)) %>%
  group_by(iso3c, country) %>%
  slice_max(date, n = 1) %>%          # último año con ambos datos
  ungroup()

# Tamaño total del comercio y balanza relativa
datos_db$log_total <- log(datos_db$Importaciones_USD + datos_db$Exportaciones_USD + 1)
datos_db$log_ratio <- log(datos_db$Exportaciones_USD + 1) - log(datos_db$Importaciones_USD + 1)

# Correlaciones para justificar la transformación
cor(log(datos_db$Exportaciones_USD + 1), log(datos_db$Importaciones_USD + 1))

## [1] 0.9943941

cor(datos_db$log_total, datos_db$log_ratio)

## [1] -0.04656584

# Estandarización Z
datos_escalados <- as.data.frame(scale(datos_db[, c("log_total", "log_ratio")]))
rownames(datos_escalados) <- datos_db$iso3c

nrow(datos_escalados)   # países que quedaron

## [1] 20

Parámetros y curva k-NN

Con 2 variables, la regla minPts ≥ 2 × m da minPts = 4.

minPts_db <- 4
epsilon_db <- 0.7

kNNdistplot(datos_escalados, k = minPts_db - 1)
abline(h = epsilon_db, col = "red", lty = 2)
title(paste("Gráfico de k-vecinos - Epsilon elegido =", epsilon_db))

set.seed(123)
resultado_db <- dbscan(datos_escalados, eps = epsilon_db, minPts = minPts_db)

datos_db$grupo <- factor(ifelse(resultado_db$cluster == 0, "Ruido",
                                paste("Grupo", resultado_db$cluster)))

table(datos_db$grupo)

## 
## Grupo 1 Grupo 2   Ruido 
##       6      11       3

Prueba de sensibilidad

valores_eps <- seq(0.2, 1.2, by = 0.1)

sensibilidad <- do.call(rbind, lapply(valores_eps, function(e) {
  r <- dbscan(datos_escalados, eps = e, minPts = minPts_db)
  data.frame(
    eps = e,
    clusteres = max(r$cluster),
    ruido = sum(r$cluster == 0),
    pct_ruido = round(100 * sum(r$cluster == 0) / nrow(datos_escalados), 1)
  )
}))

sensibilidad %>%
  kable(
    col.names = c("Epsilon", "Clústeres", "Países en ruido", "% de ruido"),
    caption = paste("Sensibilidad del resultado al valor de epsilon (minPts =", minPts_db, ")")
  ) %>%
  kable_classic(html_font = "Times New Roman", font_size = 13)

Sensibilidad del resultado al valor de epsilon (minPts = 4 )
Epsilon	Clústeres	Países en ruido	% de ruido
0.2	0	20	100
0.3	1	15	75
0.4	2	10	50
0.5	1	10	50
0.6	1	9	45
0.7	2	3	15
0.8	2	3	15
0.9	2	3	15
1.0	1	3	15
1.1	1	3	15
1.2	1	3	15

Parte 4

Grafica K-Means

fviz_cluster(
  resultado_kmeans,
  data = matriz_tipificada,
  palette = c("#0072B2", "#D55E00", "#009E73", "#CC79A7", "#E69F00", "#56B4E9"),
  ellipse.type = "norm",
  star.plot = TRUE,
  repel = TRUE,
  labelsize = 8,
  main = paste("K-Means: Clústeres de países seleccionados (k =", k_optimo, ")"),
  ggtheme = theme_minimal()
)

Grafica K-Medoides (PAM)

fviz_cluster(
  modelo_pam,
  data = matriz_tipificada,
  palette = c("#0072B2", "#D55E00", "#009E73", "#CC79A7", "#E69F00", "#56B4E9", "#555555"),
  ellipse.type = "norm",
  repel = TRUE,
  labelsize = 8,
  main = paste("K-Medoides (PAM): Clústeres de países seleccionados (k =", k_op, ")"),
  ggtheme = theme_minimal()
)

Dendrograma con el corte de Método de Ward

paleta <- c("#0072B2", "#D55E00", "#009E73", "#CC79A7", "#E69F00", "#56B4E9", "#555555")

fviz_dend(
  modelo_ward,
  k = k_ward,
  cex = 0.7,
  k_colors = paleta[1:k_ward],
  rect = TRUE,
  rect_fill = TRUE,
  main = paste0("Dendrograma - Método de Ward (k = ", k_ward, ")")
)

Grafica DBSCAN

datos_grafico <- datos_escalados
datos_grafico$iso3c <- rownames(datos_escalados)
datos_grafico$grupo <- datos_db$grupo

# Valor Z que corresponde a exportaciones = importaciones (log_ratio = 0)
equilibrio <- (0 - mean(datos_db$log_ratio)) / sd(datos_db$log_ratio)

ggplot(datos_grafico, aes(x = log_total, y = log_ratio, color = grupo)) +
  geom_point(alpha = 0.8, size = 2.8) +
  geom_text_repel(aes(label = iso3c), size = 3, show.legend = FALSE) +
  geom_hline(yintercept = equilibrio, linetype = "dashed", color = "gray50") +
  scale_color_manual(values = c("Ruido"   = "#555555",
                                "Grupo 1" = "#0072B2",
                                "Grupo 2" = "#D55E00",
                                "Grupo 3" = "#009E73",
                                "Grupo 4" = "#CC79A7")) +
  labs(
    title = "DBSCAN - Tamaño y balanza comercial de los países seleccionados",
    subtitle = paste("Parámetros: eps =", round(epsilon_db, 3), "| minPts =", minPts_db,
                     "| Línea punteada: exportaciones = importaciones"),
    x = "Tamaño del Comercio (Estandarizado Z)",
    y = "Balanza Relativa Exp-Imp (Estandarizada Z)",
    color = "Clasificación"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

Epsilon	Clústeres	Países en ruido	% de ruido
0.2	0	20	100
0.3	1	15	75
0.4	2	10	50
0.5	1	10	50
0.6	1	9	45
0.7	2	3	15
0.8	2	3	15
0.9	2	3	15
1.0	1	3	15
1.1	1	3	15
1.2	1	3	15

Epsilon	Clústeres	Países en ruido	% de ruido
0.2	0	20	100
0.3	1	15	75
0.4	2	10	50
0.5	1	10	50
0.6	1	9	45
0.7	2	3	15
0.8	2	3	15
0.9	2	3	15
1.0	1	3	15
1.1	1	3	15
1.2	1	3	15

Epsilon	Clústeres	Países en ruido	% de ruido
0.2	0	20	100
0.3	1	15	75
0.4	2	10	50
0.5	1	10	50
0.6	1	9	45
0.7	2	3	15
0.8	2	3	15
0.9	2	3	15
1.0	1	3	15
1.1	1	3	15
1.2	1	3	15