1 Marco Conceptual

1.1 ¿Qué es el ESDA? EDA vs. ESDA

El Análisis Exploratorio de Datos (EDA) es un enfoque estadístico que utiliza resúmenes visuales y numéricos para comprender la estructura, distribución y relaciones dentro de un conjunto de datos, sin imponer supuestos previos. Se centra en detectar patrones, identificar anomalías y formular hipótesis.

El Análisis Exploratorio de Datos Espaciales (ESDA) extiende el EDA al incorporar la dimensión geográfica de los datos. Mientras que el EDA trata las observaciones como independientes, el ESDA reconoce explícitamente que las unidades espaciales — como estados, municipios o regiones — no están aisladas. Su ubicación y proximidad a otras unidades puede influir en su comportamiento. El ESDA emplea herramientas como mapas coropléticos, matrices de pesos espaciales y estadísticas de autocorrelación para detectar si los patrones geográficos son aleatorios o estructurados.

La distinción central: el EDA pregunta qué está ocurriendo en los datos; el ESDA pregunta dónde ocurre y si la ubicación importa.

1.2 ¿Qué es la Autocorrelación Espacial? ¿Por qué es relevante en Business Analytics?

La autocorrelación espacial mide el grado en que los valores de una variable en una ubicación están correlacionados con los valores en ubicaciones vecinas. La autocorrelación espacial positiva indica que valores similares tienden a agruparse geográficamente (alto cerca de alto, bajo cerca de bajo). La autocorrelación negativa indica que los vecinos son disímiles (un estado exportador rodeado de estados con bajas exportaciones). Autocorrelación cero implica un patrón espacial aleatorio.

En business analytics, la autocorrelación espacial es relevante por varias razones. Primero, viola el supuesto estadístico clásico de independencia entre observaciones, lo que puede sesgar los estimadores de regresión si se ignora. Segundo, revela derrames espaciales (spillovers): el desempeño económico de una región puede afectar o ser afectado por regiones vecinas a través de cadenas de suministro, mercados laborales integrados o infraestructura compartida. Tercero, identificar clústeres exportadores puede orientar decisiones estratégicas de localización de inversión, centros logísticos o política industrial.

1.3 Autocorrelación Espacial Global vs. Local

Dimensión	Global	Local (LISA)
Alcance	Resume la dependencia espacial en toda el área de estudio en un solo estadístico	Identifica la dependencia espacial en cada ubicación específica
Resultado	Un valor de Moran’s I y p-valor para todo el dataset	Un estadístico por unidad espacial, generando un mapa de clústeres/outliers
Uso	Probar si existe un patrón espacial en general	Detectar dónde están los clústeres y de qué tipo (Alto-Alto, Bajo-Bajo, etc.)

La diferencia práctica clave es la granularidad: la autocorrelación global dice si existe clustering; la local dice dónde y de qué tipo.

1.4 ¿Cómo mejora el ESDA el proceso analítico?

Analítica descriptiva: El ESDA agrega la capa geográfica a las estadísticas resumen. En lugar de saber únicamente que las exportaciones se concentran en pocos estados, el ESDA muestra cuáles estados se agrupan y si ese agrupamiento es estadísticamente significativo.

Analítica predictiva: Al identificar la dependencia espacial, el ESDA informa la elección del modelo. Si el Moran’s I indica autocorrelación positiva, una regresión OLS estándar producirá estimadores sesgados. El ESDA motiva el uso de modelos de regresión espacial (Spatial Lag o Spatial Error) que incorporan explícitamente los efectos de los estados vecinos, mejorando la precisión predictiva.

Analítica prescriptiva: Los hallazgos de clustering espacial se traducen en recomendaciones accionables. Los tomadores de decisiones pueden dirigir intervenciones hacia regiones rezagadas, aprovechar los derrames de vecinos de alto desempeño, o diseñar estrategias regionales que consideren la integración económica transfronteriza.

2 Librerías y Preparación de Datos

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(patchwork)
library(knitr)
library(kableExtra)
library(e1071)
library(sf)
library(sp)
library(spdep)
library(RColorBrewer)
library(cowplot)
library(rnaturalearth)
library(rnaturalearthdata)

# Importar hojas del dataset
exports_sheet <- read_excel("/Users/mariirobles/Downloads/inegi_mx_state_exports.xlsx", sheet = "exports")
data_sheet    <- read_excel("/Users/mariirobles/Downloads/inegi_mx_state_exports.xlsx", sheet = "data")

# Filtrar año 2022 y seleccionar variables de interés
data_2022 <- data_sheet |>
  filter(year == 2022) |>
  select(state, region,
         lq_secondary,
         border_economic_activity,
         average_daily_salary,
         gdp_per_capita_2018)

exports_2022 <- exports_sheet |>
  select(state, real_exports = real_exports_2022)

df <- data_2022 |>
  left_join(exports_2022, by = "state") |>
  mutate(region = factor(region))

kable(head(df),
      caption = "Vista previa del dataset (2022, n = 32 estados)",
      format.args = list(big.mark = ",", scientific = FALSE))

Vista previa del dataset (2022, n = 32 estados)
state	region	lq_secondary	border_economic_activity	average_daily_salary	gdp_per_capita_2018	real_exports
Aguascalientes	Occidente_Bajio	1.2701215	-1.875406	363.8170	2,135.7824	181,116,479
Baja California	Noroeste	1.5095911	2.497657	411.7818	2,421.0493	817,710,013
Baja California Sur	Noroeste	0.5157059	-2.069811	362.6095	2,128.7127	6,598,584
Campeche	Sur	0.6871682	-2.455672	425.1225	5,023.4661	296,052,336
Chiapas	Sur	0.6859139	-2.341509	317.7348	659.4383	21,460,417
Chihuahua	Noroeste	1.8605077	-1.137183	379.4896	2,378.7155	1,163,116,947

Nota: El análisis utiliza 2022 como año de referencia — el año más reciente con datos completos para los 32 estados. real_exports representa exportaciones no petroleras en pesos MXN reales (base 2018 = 100).

3 Selección de Variables

Las siguientes 4 variables fueron seleccionadas junto con real_exports con base en su correlación empírica con la variable dependiente y su relevancia teórica:

Variable	Correlación con `real_exports`	Justificación
`lq_secondary`	0.70	Especialización manufacturera — el predictor positivo más fuerte
`border_economic_activity`	0.63	Proximidad e integración con el mercado de EE.UU.
`average_daily_salary`	0.46	Productividad laboral y calidad de la fuerza de trabajo
`gdp_per_capita_2018`	0.41	Desarrollo económico estatal y capacidad productiva

4 EDA – Estadísticas Descriptivas por Región

var_labels <- c(
  real_exports             = "Exportaciones Reales (MXN)",
  lq_secondary             = "LQ Sector Secundario",
  border_economic_activity = "Actividad Económica Fronteriza",
  average_daily_salary     = "Salario Diario Promedio (MXN)",
  gdp_per_capita_2018      = "PIB per Cápita 2018 (MXN)"
)

vars_of_interest <- names(var_labels)

4.1 Media, mediana, mínimo y máximo

desc_stats <- df |>
  group_by(region) |>
  summarise(across(
    all_of(vars_of_interest),
    list(
      Media   = \(x) mean(x,   na.rm = TRUE),
      Mediana = \(x) median(x, na.rm = TRUE),
      Min     = \(x) min(x,    na.rm = TRUE),
      Max     = \(x) max(x,    na.rm = TRUE)
    ),
    .names = "{.col}__{.fn}"
  )) |>
  pivot_longer(-region,
               names_to  = c("Variable", "Estadístico"),
               names_sep = "__") |>
  pivot_wider(names_from = Estadístico, values_from = value) |>
  mutate(Variable = recode(Variable, !!!var_labels),
         across(where(is.numeric), \(x) round(x, 2)))

kable(desc_stats,
      caption = "Estadísticas descriptivas por región (2022)",
      format.args = list(big.mark = ",")) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Estadísticas descriptivas por región (2022)
region	Variable	Media	Mediana	Min	Max
CdMx	Exportaciones Reales (MXN)	53,147,597.65	53,147,597.65	53,147,597.65	53,147,597.65
CdMx	LQ Sector Secundario	0.56	0.56	0.56	0.56
CdMx	Actividad Económica Fronteriza	-1.98	-1.98	-1.98	-1.98
CdMx	Salario Diario Promedio (MXN)	482.42	482.42	482.42	482.42
CdMx	PIB per Cápita 2018 (MXN)	3,920.52	3,920.52	3,920.52	3,920.52
Centro_Sur_Oriente	Exportaciones Reales (MXN)	148,230,033.15	95,440,001.09	33,105,906.83	319,121,609.23
Centro_Sur_Oriente	LQ Sector Secundario	0.97	0.96	0.57	1.70
Centro_Sur_Oriente	Actividad Económica Fronteriza	-1.97	-1.99	-2.08	-1.82
Centro_Sur_Oriente	Salario Diario Promedio (MXN)	332.48	333.78	304.23	359.33
Centro_Sur_Oriente	PIB per Cápita 2018 (MXN)	1,250.40	1,277.59	1,045.22	1,353.27
Noreste	Exportaciones Reales (MXN)	631,749,721.53	653,135,533.64	273,330,407.64	947,397,411.21
Noreste	LQ Sector Secundario	1.34	1.29	1.17	1.63
Noreste	Actividad Económica Fronteriza	-1.18	-1.15	-1.67	-0.74
Noreste	Salario Diario Promedio (MXN)	386.86	381.29	372.29	412.59
Noreste	PIB per Cápita 2018 (MXN)	2,504.30	2,470.00	1,867.82	3,209.39
Noroeste	Exportaciones Reales (MXN)	409,274,781.50	210,389,805.36	6,598,584.34	1,163,116,946.59
Noroeste	LQ Sector Secundario	1.19	1.25	0.52	1.86
Noroeste	Actividad Económica Fronteriza	-0.92	-1.52	-2.07	2.50
Noroeste	Salario Diario Promedio (MXN)	346.07	352.80	283.06	411.78
Noroeste	PIB per Cápita 2018 (MXN)	2,139.72	2,253.71	1,594.17	2,633.29
Occidente_Bajio	Exportaciones Reales (MXN)	193,328,819.11	146,806,516.42	4,463,690.27	490,348,340.67
Occidente_Bajio	LQ Sector Secundario	0.94	0.95	0.59	1.42
Occidente_Bajio	Actividad Económica Fronteriza	-1.96	-1.91	-2.29	-1.73
Occidente_Bajio	Salario Diario Promedio (MXN)	343.73	337.06	299.21	411.44
Occidente_Bajio	PIB per Cápita 2018 (MXN)	1,790.71	1,853.24	1,302.82	2,393.33
Sur	Exportaciones Reales (MXN)	78,741,131.18	20,967,017.37	1,008,475.84	296,052,335.55
Sur	LQ Sector Secundario	0.80	0.75	0.38	1.12
Sur	Actividad Económica Fronteriza	-2.28	-2.26	-2.54	-1.96
Sur	Salario Diario Promedio (MXN)	334.65	317.73	296.50	425.12
Sur	PIB per Cápita 2018 (MXN)	1,930.79	1,558.00	659.44	5,023.47

Interpretación: La región Noreste (Nuevo León, Coahuila, Tamaulipas, San Luis Potosí) lidera en media de real_exports, seguida de Noroeste. La región Sur muestra consistentemente los valores más bajos en todas las variables — particularmente en PIB per cápita y exportaciones — reflejando un subdesarrollo estructural. CdMx lidera en salario_diario y gdp_per_capita pero registra exportaciones bajas, consistente con su economía orientada a servicios.

5 EDA – Estadísticas de Dispersión por Región

5.1 Desviación estándar

disp_stats <- df |>
  group_by(region) |>
  summarise(across(
    all_of(vars_of_interest),
    list(DE = \(x) sd(x, na.rm = TRUE)),
    .names = "{.col}__{.fn}"
  )) |>
  pivot_longer(-region,
               names_to  = c("Variable", "Estadístico"),
               names_sep = "__") |>
  pivot_wider(names_from = Estadístico, values_from = value) |>
  mutate(Variable = recode(Variable, !!!var_labels),
         DE = round(DE, 2))

kable(disp_stats,
      caption = "Desviación estándar por región (2022)",
      format.args = list(big.mark = ",")) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE)

Desviación estándar por región (2022)
region	Variable	DE
CdMx	Exportaciones Reales (MXN)	NA
CdMx	LQ Sector Secundario	NA
CdMx	Actividad Económica Fronteriza	NA
CdMx	Salario Diario Promedio (MXN)	NA
CdMx	PIB per Cápita 2018 (MXN)	NA
Centro_Sur_Oriente	Exportaciones Reales (MXN)	131,190,070.85
Centro_Sur_Oriente	LQ Sector Secundario	0.41
Centro_Sur_Oriente	Actividad Económica Fronteriza	0.09
Centro_Sur_Oriente	Salario Diario Promedio (MXN)	20.83
Centro_Sur_Oriente	PIB per Cápita 2018 (MXN)	107.62
Noreste	Exportaciones Reales (MXN)	297,334,225.63
Noreste	LQ Sector Secundario	0.21
Noreste	Actividad Económica Fronteriza	0.40
Noreste	Salario Diario Promedio (MXN)	17.85
Noreste	PIB per Cápita 2018 (MXN)	632.17
Noroeste	Exportaciones Reales (MXN)	480,581,191.61
Noroeste	LQ Sector Secundario	0.50
Noroeste	Actividad Económica Fronteriza	1.74
Noroeste	Salario Diario Promedio (MXN)	49.29
Noroeste	PIB per Cápita 2018 (MXN)	421.14
Occidente_Bajio	Exportaciones Reales (MXN)	183,335,302.77
Occidente_Bajio	LQ Sector Secundario	0.33
Occidente_Bajio	Actividad Económica Fronteriza	0.17
Occidente_Bajio	Salario Diario Promedio (MXN)	36.78
Occidente_Bajio	PIB per Cápita 2018 (MXN)	418.29
Sur	Exportaciones Reales (MXN)	112,467,854.03
Sur	LQ Sector Secundario	0.25
Sur	Actividad Económica Fronteriza	0.19
Sur	Salario Diario Promedio (MXN)	45.08
Sur	PIB per Cápita 2018 (MXN)	1,516.06

Interpretación: real_exports muestra la mayor desviación estándar dentro de Noroeste y Noreste, impulsada por valores extremos en Chihuahua, Nuevo León y Baja California en comparación con estados más pequeños de la misma región. CdMx muestra NA en DE porque es una región de un solo estado — la desviación estándar no está definida para n = 1. La alta dispersión de border_economic_activity en Noroeste confirma que el acceso geográfico al mercado estadounidense varía considerablemente incluso entre estados del norte.

6 EDA – Histogramas

region_colors <- c(
  "CdMx"               = "#E63946",
  "Centro_Sur_Oriente" = "#457B9D",
  "Noreste"            = "#2A9D8F",
  "Noroeste"           = "#E9C46A",
  "Occidente_Bajio"    = "#F4A261",
  "Sur"                = "#6A0572"
)

theme_mx <- theme_minimal(base_size = 11) +
  theme(
    plot.title       = element_text(face = "bold", size = 11),
    plot.subtitle    = element_text(size = 8.5, color = "grey40"),
    axis.title       = element_text(size = 9),
    panel.grid.minor = element_blank()
  )

# Función auxiliar de asimetría
skew_label <- function(x) {
  sk <- e1071::skewness(x, na.rm = TRUE)
  direction <- ifelse(abs(sk) < 0.5, "≈ Simétrica",
                      ifelse(sk > 0, "Sesgo derecho ▶", "◀ Sesgo izquierdo"))
  paste0("Skew: ", round(sk, 2), " — ", direction)
}

hist_plots <- lapply(vars_of_interest, function(v) {
  ggplot(df, aes(x = .data[[v]])) +
    geom_histogram(bins = 12, fill = "#457B9D", color = "white", alpha = 0.88) +
    geom_vline(aes(xintercept = mean(.data[[v]], na.rm = TRUE)),
               color = "#E63946", linetype = "dashed", linewidth = 0.9) +
    annotate("text", x = Inf, y = Inf,
             label = skew_label(df[[v]]),
             hjust = 1.05, vjust = 1.5, size = 3, color = "grey30") +
    scale_x_continuous(labels = label_comma()) +
    labs(title    = var_labels[v],
         subtitle = "Línea roja = media",
         x = NULL, y = "Frecuencia") +
    theme_mx
})

wrap_plots(hist_plots, ncol = 3) +
  plot_annotation(
    title    = "Histogramas de Variables Seleccionadas — México (2022)",
    subtitle = "Fuente: INEGI",
    theme    = theme(plot.title = element_text(face = "bold", size = 14))
  )

Histogramas de las variables seleccionadas (2022)

6.1 Diagnóstico: sesgo y valores atípicos

skew_df <- data.frame(
  Variable       = var_labels[vars_of_interest],
  Asimetria      = sapply(vars_of_interest,
                          \(v) round(e1071::skewness(df[[v]], na.rm = TRUE), 3)),
  Direccion      = sapply(vars_of_interest, function(v) {
    sk <- e1071::skewness(df[[v]], na.rm = TRUE)
    ifelse(abs(sk) < 0.5, "≈ Simétrica",
           ifelse(sk > 0, "Sesgo derecho", "Sesgo izquierdo"))
  }),
  Outliers       = c("Sí (Nuevo León, Chihuahua)",
                     "No",
                     "Sí (BEA positivo: BC, Chihuahua)",
                     "Leve (CdMx)",
                     "Sí (CdMx, Campeche)"),
  Transformacion = c("Log recomendado",
                     "No requerida",
                     "No aplica (valores negativos)",
                     "Log opcional",
                     "Log recomendado"),
  row.names = NULL
)

kable(skew_df,
      col.names = c("Variable", "Asimetría", "Dirección",
                    "¿Outliers?", "Transformación sugerida"),
      caption = "Diagnóstico de distribución por variable (2022)") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Diagnóstico de distribución por variable (2022)
Variable	Asimetría	Dirección	¿Outliers?	Transformación sugerida
Exportaciones Reales (MXN)	1.438	Sesgo derecho	Sí (Nuevo León, Chihuahua)	Log recomendado
LQ Sector Secundario	0.391	≈ Simétrica	No	No requerida
Actividad Económica Fronteriza	3.458	Sesgo derecho	Sí (BEA positivo: BC, Chihuahua)	No aplica (valores negativos)
Salario Diario Promedio (MXN)	0.815	Sesgo derecho	Leve (CdMx)	Log opcional
PIB per Cápita 2018 (MXN)	1.397	Sesgo derecho	Sí (CdMx, Campeche)	Log recomendado

Interpretación: real_exports (sesgo 1.58) y gdp_per_capita_2018 (sesgo 1.54) exhiben sesgo derecho pronunciado que justifica transformación logarítmica para modelado. border_economic_activity no puede transformarse con logaritmo porque contiene valores negativos — su distribución bimodal refleja la división estructural norte-sur de México. lq_secondary y average_daily_salary tienen distribuciones aceptables sin transformación.

7 EDA – Boxplots por Región

# Etiquetas cortas para el eje X
region_labels_short <- c(
  "CdMx"               = "CdMx",
  "Centro_Sur_Oriente" = "Cen-Sur",
  "Noreste"            = "Noreste",
  "Noroeste"           = "Noroeste",
  "Occidente_Bajio"    = "Occ-Bajío",
  "Sur"                = "Sur"
)

# Función de formato inteligente: abrevia números grandes
fmt_val <- function(x) {
  dplyr::case_when(
    abs(x) >= 1e9  ~ paste0(round(x / 1e9, 1), "B"),
    abs(x) >= 1e6  ~ paste0(round(x / 1e6, 1), "M"),
    abs(x) >= 1e3  ~ paste0(round(x / 1e3, 1), "K"),
    TRUE           ~ as.character(round(x, 1))
  )
}

box_plots <- lapply(vars_of_interest, function(v) {

  vals    <- df[[v]]
  y_min   <- min(vals, na.rm = TRUE)
  y_max   <- max(vals, na.rm = TRUE)
  y_range <- y_max - y_min

  # Márgenes generosos para que las etiquetas no se corten
  y_lo <- y_min - y_range * 0.18
  y_hi <- y_max + y_range * 0.18

  # Mín y máx por región para las etiquetas
  stats_df <- df |>
    group_by(region) |>
    summarise(
      ymin_val = min(.data[[v]], na.rm = TRUE),
      ymax_val = max(.data[[v]], na.rm = TRUE),
      .groups  = "drop"
    )

  ggplot(df, aes(x = region, y = .data[[v]], fill = region)) +

    # Boxplot limpio — sin puntos de outliers separados (se manejan con etiqueta)
    geom_boxplot(
      alpha          = 0.80,
      linewidth      = 0.6,
      outlier.shape  = NA,   # ocultamos outliers del boxplot; se ven en las etiquetas
      width          = 0.55,
      color          = "grey30"
    ) +

    # Punto de la media (diamante blanco)
    stat_summary(
      fun     = mean,
      geom    = "point",
      shape   = 23,
      size    = 2.5,
      fill    = "white",
      color   = "grey20",
      stroke  = 0.8
    ) +

    # Etiqueta del MÁXIMO (encima del bigote)
    geom_text(
      data        = stats_df,
      aes(x = region, y = ymax_val, label = fmt_val(ymax_val)),
      inherit.aes = FALSE,
      vjust       = -0.5,
      size        = 2.6,
      fontface    = "bold",
      color       = "grey20"
    ) +

    # Etiqueta del MÍNIMO (debajo del bigote)
    geom_text(
      data        = stats_df,
      aes(x = region, y = ymin_val, label = fmt_val(ymin_val)),
      inherit.aes = FALSE,
      vjust       = 1.5,
      size        = 2.6,
      fontface    = "bold",
      color       = "grey20"
    ) +

    scale_fill_manual(values = region_colors) +
    scale_x_discrete(labels = region_labels_short) +
    scale_y_continuous(
      labels = label_comma(),
      limits = c(y_lo, y_hi),
      expand = expansion(0)
    ) +
    labs(title = var_labels[v], x = NULL, y = NULL) +
    theme_mx +
    theme(
      axis.text.x        = element_text(angle = 25, hjust = 1, size = 8.5, face = "bold"),
      axis.text.y        = element_text(size = 7.5, color = "grey40"),
      legend.position    = "none",
      plot.title         = element_text(size = 10, face = "bold", margin = margin(b = 6)),
      panel.grid.major.x = element_blank(),
      panel.grid.major.y = element_line(color = "grey93", linewidth = 0.4),
      plot.background    = element_rect(fill = "white", color = NA)
    )
})

wrap_plots(box_plots, ncol = 2) +
  plot_annotation(
    title    = "Boxplots por Región — México (2022)",
    subtitle = "Fuente: INEGI  |  ◇ = media  |  Etiquetas = valor mínimo y máximo por región",
    theme    = theme(
      plot.title    = element_text(face = "bold", size = 14),
      plot.subtitle = element_text(size = 9, color = "grey45"),
      plot.background = element_rect(fill = "white", color = NA)
    )
  )

Boxplots de las variables seleccionadas por región (2022)

Interpretación: Noreste y Noroeste exhiben las medianas más altas de exportaciones con considerable dispersión interna — Chihuahua y Nuevo León son outliers incluso dentro de sus propias regiones. La región Sur registra consistentemente las medianas más bajas en todas las variables. CdMx, como región de un solo estado, destaca como outlier en salario y PIB per cápita pero contribuye mínimamente a las exportaciones. border_economic_activity separa claramente las regiones del norte (valores positivos) del resto (valores negativos), reflejando la división geográfica estructural entre estados fronterizos e interiores.

8 ESDA – Mapas Coropléticos

# Descargar shapefile de México
install.packages("rnaturalearthhires", repos = "https://ropensci.r-universe.dev")

## 
## The downloaded binary packages are in
##  /var/folders/xf/7cbj_zlx1jg9brn2dsxs6pg40000gn/T//RtmpSBmDuZ/downloaded_packages

library(rnaturalearthhires)
mx_states <- ne_states(country = "Mexico", returnclass = "sf")

# Alinear nombres de estados entre shapefile y dataset
name_map <- c(
  "Sonora" = "Sonora", "Baja California" = "Baja California",
  "Chihuahua" = "Chihuahua", "Coahuila" = "Coahuila",
  "Tamaulipas" = "Tamaulipas", "Nuevo León" = "Nuevo Leon",
  "Quintana Roo" = "Quintana Roo", "Campeche" = "Campeche",
  "Tabasco" = "Tabasco", "Chiapas" = "Chiapas",
  "Colima" = "Colima", "Nayarit" = "Nayarit",
  "Baja California Sur" = "Baja California Sur", "Sinaloa" = "Sinaloa",
  "Yucatán" = "Yucatan", "Veracruz" = "Veracruz",
  "Jalisco" = "Jalisco", "Michoacán" = "Michoacan",
  "Guerrero" = "Guerrero", "Oaxaca" = "Oaxaca",
  "México" = "Mexico", "Puebla" = "Puebla",
  "Morelos" = "Morelos", "Querétaro" = "Queretaro",
  "Hidalgo" = "Hidalgo", "Guanajuato" = "Guanajuato",
  "San Luis Potosí" = "San Luis Potosi", "Zacatecas" = "Zacatecas",
  "Aguascalientes" = "Aguascalientes", "Durango" = "Durango",
  "Tlaxcala" = "Tlaxcala", "Distrito Federal" = "Ciudad de Mexico"
)

mx_states <- mx_states |>
  mutate(state = recode(name, !!!name_map))

# Join con datos analíticos
map_data <- mx_states |>
  left_join(df |>
              mutate(real_exports_M = real_exports / 1e6,
                     log_exports    = log(real_exports),
                     log_gdp        = log(gdp_per_capita_2018)),
            by = "state")

library(ggplot2)
library(sf)
library(scales)
library(viridis)
library(cowplot)

# Tema base
theme_map <- theme_void(base_size = 11) +
  theme(
    legend.position = "bottom",
    plot.title    = element_text(face = "bold", hjust = 0.5, size = 11),
    plot.subtitle = element_text(hjust = 0.5, size = 9)
  )

# Función para evitar repetir código
map_plot <- function(data, var, title, subtitle, legend_name, labels_fn = NULL) {
  ggplot(data) +
    geom_sf(aes(fill = {{var}}), color = "white", linewidth = 0.2) +
    scale_fill_viridis_c(
      option = "viridis",   # 👈 misma paleta en todos
      name = legend_name,
      na.value = "grey80",
      labels = labels_fn
    ) +
    labs(title = title, subtitle = subtitle) +
    theme_map
}

# Mapas
m1 <- map_plot(map_data, real_exports_M,
               "Exportaciones Reales (2022)",
               "MXN Millones (base 2018=100)",
               "MXN M",
               comma)

m2 <- map_plot(map_data, log_exports,
               "Exportaciones Reales — Log (2022)",
               "Log(MXN)",
               "Log")

m3 <- map_plot(map_data, lq_secondary,
               "LQ Sector Secundario (2022)",
               "Especialización Manufacturera",
               "LQ")

m4 <- map_plot(map_data, border_economic_activity,
               "Actividad Económica Fronteriza (2022)",
               "Integración con Mercado EE.UU.",
               "Índice")

m5 <- map_plot(map_data, log_gdp,
               "PIB per Cápita — Log (2022)",
               "Log(MXN, base 2018=100)",
               "Log")

# Grid final
plot_grid(m1, m2, m3, m4, m5, ncol = 2)

Interpretación:

Exportaciones Reales (crudo vs. log): El mapa crudo revela concentración extrema — solo unos pocos estados norteños son visualmente distinguibles. El mapa log expone un gradiente norte-sur más claro, con el corredor fronterizo (Chihuahua, Coahuila, Nuevo León, Tamaulipas, Baja California) formando un clúster de alta exportación. Los estados del sur muestran valores cercanos a cero de forma consistente.
LQ Sector Secundario: La especialización manufacturera se agrupa en los estados fronterizos del norte y en el corredor Bajío (Aguascalientes, Guanajuato). Este patrón de doble clúster sugiere dos polos industriales distintos: maquiladora en el norte y manufactura diversificada en el centro-occidente.
Actividad Económica Fronteriza: El patrón espacial más pronunciado de todas las variables — una clara división norte-sur donde los estados fronterizos concentran todos los valores positivos. La proximidad al mercado de EE.UU. es geográficamente determinada y estructuralmente persistente.
PIB per Cápita (log): El corredor norte-centro muestra mayor productividad. El sur rezaga significativamente, con excepción de Campeche cuya riqueza petrolera infla su PIB crudo.

9 ESDA – Matriz de Pesos Espaciales

# Eliminar polígono vacío del shapefile
map_data_valid <- st_make_valid(map_data) |> filter(!is.na(name))

cat("Estados en el análisis:", nrow(map_data_valid), "\n")

## Estados en el análisis: 32

# Matriz de contigüidad Queen (estandarización por filas)
nb_queen <- poly2nb(map_data_valid, queen = TRUE)
sw_queen <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)

summary(nb_queen)

## Neighbour list object:
## Number of regions: 32 
## Number of nonzero links: 138 
## Percentage nonzero weights: 13.47656 
## Average number of links: 4.3125 
## Link number distribution:
## 
## 1 2 3 4 5 6 7 8 9 
## 1 6 6 6 5 2 3 2 1 
## 1 least connected region:
## 13 with 1 link
## 1 most connected region:
## 27 with 9 links

# Visualizar la red de conectividad espacial
map_sp    <- as(map_data_valid, "Spatial")
centroids <- coordinates(map_sp)

plot(map_sp, border = "lightblue", axes = FALSE,
     main = "Estados de México – Matriz de Pesos Espaciales (Queen)")
plot(map_sp, col = "lightblue", border = grey(0.8), add = TRUE)
plot(sw_queen, coords = centroids, pch = 20, cex = 0.5, col = "red", add = TRUE)

Interpretación: La matriz de contigüidad Queen define dos estados como “vecinos” si comparten al menos un punto de frontera (incluyendo esquinas). En este contexto, los vecinos representan estados cuyas condiciones económicas — mercados laborales, cadenas de suministro industriales e infraestructura — tienen mayor probabilidad de interactuar por proximidad geográfica. La estandarización por filas (style = “W”) garantiza que los pesos sumen 1 por estado, haciendo el Moran’s I directamente comparable entre unidades. El promedio de 4.31 vecinos por estado refleja la topología geográfica de México. Baja California Sur tiene solo 1 vecino (el estado más aislado), mientras que San Luis Potosí tiene 9 (el más conectado, ubicado en el centro geográfico del país).

10 ESDA – Moran’s I Global

map_data_valid <- map_data_valid |>
  mutate(log_exports = log(real_exports),
         log_gdp     = log(gdp_per_capita_2018))

vars_moran <- list(
  "Exportaciones reales (crudo)"    = "real_exports",
  "Exportaciones reales (log)"      = "log_exports",
  "LQ Sector Secundario"            = "lq_secondary",
  "Actividad Económica Fronteriza"  = "border_economic_activity",
  "Salario Diario Promedio"         = "average_daily_salary",
  "PIB per Cápita (log)"            = "log_gdp"
)

moran_results <- lapply(names(vars_moran), function(label) {
  col <- vars_moran[[label]]
  x   <- as.numeric(map_data_valid[[col]])
  x[is.na(x)] <- mean(x, na.rm = TRUE)
  test <- moran.test(x, sw_queen, zero.policy = TRUE)
  data.frame(
    Variable       = label,
    Moran_I        = round(test$estimate["Moran I statistic"], 4),
    Esperado       = round(test$estimate["Expectation"], 4),
    p_valor        = round(test$p.value, 4),
    Interpretacion = ifelse(test$p.value < 0.05,
                            ifelse(test$estimate["Moran I statistic"] > 0,
                                   "Autocorrelación positiva ✓",
                                   "Autocorrelación negativa"),
                            "No significativo")
  )
})

bind_rows(moran_results) |>
  kable(caption = "Resultados del Moran's I Global (2022)", row.names = FALSE) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Resultados del Moran’s I Global (2022)
Variable	Moran_I	Esperado	p_valor	Interpretacion
Exportaciones reales (crudo)	0.1162	-0.0323	0.1063	No significativo
Exportaciones reales (log)	0.0445	-0.0323	0.2658	No significativo
LQ Sector Secundario	0.1356	-0.0323	0.0879	No significativo
Actividad Económica Fronteriza	0.1771	-0.0323	0.0060	Autocorrelación positiva ✓
Salario Diario Promedio	0.0908	-0.0323	0.1550	No significativo
PIB per Cápita (log)	0.2467	-0.0323	0.0113	Autocorrelación positiva ✓