Autocorrelación Espacial

Actividad 02

Planeación estratégica basada en analítica prescriptiva

Grupo 503

Profesor Rodolfo Miguel Gameros

Equipo 7:

A00833113 - Avril Lobato

A01771127 - Lesly Darian Romero Vázquez

A00831105 - Jazmín Cortez Mendoza

A01284611 - Lisset Hernández

Análisis Exploratorio de Datos Espaciales (ESDA: Exploratory Spatial Data Analysis)

El Análisis exploratorio de datos espacial involucra diferentes herramientas para visualizar los patrones espaciales existentes en la base de datos. Dichos datos se pueden segregar o descomponer los patrones espaciales en diferentes elementos a lo largo des mapa, y a su vez medir la extensión de estas autocorrelaciones espaciales. El principal paquete de R que se utiliza para este propósito es ‘geostan’.

Importación de bases de datos y geocercas

A continuación, se realiza a importación de la bd a nivel estatal

mx_state <- read.csv("state_data.csv")
head(n=10, mx_state)

##                  state year state_id crime_rate unemployment employment
## 1       Aguascalientes 2021     1057       6.75         0.04       0.97
## 2      Baja California 2021     2304      84.67         0.01       0.98
## 3  Baja California Sur 2021     2327       8.52         0.03       0.97
## 4             Campeche 2021     1086       9.22         0.02       0.98
## 5              Chiapas 2021     1182      10.01         0.05       0.97
## 6            Chihuahua 2021      888      71.98         0.04       0.97
## 7     Ciudad de Mexico 2021     1114      11.58         0.06       0.94
## 8             Coahuila 2021      933       4.92         0.04       0.94
## 9               Colima 2021     1111      71.97         0.05       0.97
## 10             Durango 2021     1004       7.91         0.05       0.97
##    business_activity real_wage real_ave_month_income pop_density lq_primary
## 1              -1.90    361.02               5641.67      261.21       0.16
## 2               2.47    388.22               7599.62       53.19       0.47
## 3              -2.12    345.57               8660.90       11.27       0.73
## 4              -2.44    414.48               5357.29       16.42       0.73
## 5              -2.41    312.37               6581.28       77.16       1.56
## 6              -1.25    362.93               5708.75       15.30       0.64
## 7              -2.08    466.80               4528.42     6195.59       0.03
## 8              -1.38    365.18               6480.66       21.44       0.17
## 9              -2.31    324.82               6097.68      136.06       1.19
## 10             -1.93    289.78               4969.61       15.07       1.05
##    lq_secondary lq_tertiary new_fdi_real_mxn log_new_fdi_real_mxn  region_a
## 1          1.24        1.00          1270.03                 3.10     Bajio
## 2          1.62        0.86         20407.81                 4.31     Norte
## 3          0.51        1.13         19022.55                 4.28 Occidente
## 4          0.78        1.06          2390.64                 3.38       Sur
## 5          0.88        0.99           189.58                 2.28       Sur
## 6          1.97        0.80          7153.63                 3.86     Norte
## 7          0.58        1.14         34394.95                 4.54    Centro
## 8          1.55        0.93          6254.19                 3.80     Norte
## 9          0.58        1.07            63.41                 1.80 Occidente
## 10         1.15        0.96          -483.58                -2.68 Occidente
##    region_b
## 1         2
## 2         3
## 3         4
## 4         5
## 5         5
## 6         3
## 7         1
## 8         3
## 9         4
## 10        4

glimpse(mx_state)

## Rows: 32
## Columns: 17
## $ state                 <chr> "Aguascalientes", "Baja California", "Baja Calif…
## $ year                  <int> 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, …
## $ state_id              <int> 1057, 2304, 2327, 1086, 1182, 888, 1114, 933, 11…
## $ crime_rate            <dbl> 6.75, 84.67, 8.52, 9.22, 10.01, 71.98, 11.58, 4.…
## $ unemployment          <dbl> 0.04, 0.01, 0.03, 0.02, 0.05, 0.04, 0.06, 0.04, …
## $ employment            <dbl> 0.97, 0.98, 0.97, 0.98, 0.97, 0.97, 0.94, 0.94, …
## $ business_activity     <dbl> -1.90, 2.47, -2.12, -2.44, -2.41, -1.25, -2.08, …
## $ real_wage             <dbl> 361.02, 388.22, 345.57, 414.48, 312.37, 362.93, …
## $ real_ave_month_income <dbl> 5641.67, 7599.62, 8660.90, 5357.29, 6581.28, 570…
## $ pop_density           <dbl> 261.21, 53.19, 11.27, 16.42, 77.16, 15.30, 6195.…
## $ lq_primary            <dbl> 0.16, 0.47, 0.73, 0.73, 1.56, 0.64, 0.03, 0.17, …
## $ lq_secondary          <dbl> 1.24, 1.62, 0.51, 0.78, 0.88, 1.97, 0.58, 1.55, …
## $ lq_tertiary           <dbl> 1.00, 0.86, 1.13, 1.06, 0.99, 0.80, 1.14, 0.93, …
## $ new_fdi_real_mxn      <dbl> 1270.03, 20407.81, 19022.55, 2390.64, 189.58, 71…
## $ log_new_fdi_real_mxn  <dbl> 3.10, 4.31, 4.28, 3.38, 2.28, 3.86, 4.54, 3.80, …
## $ region_a              <chr> "Bajio", "Norte", "Occidente", "Sur", "Sur", "No…
## $ region_b              <int> 2, 3, 4, 5, 5, 3, 1, 3, 4, 4, 2, 5, 1, 4, 1, 4, …

Seguidamente, se leen los archivos de las formas geométricas (geocercas) que determinarán los mapas a realizar después para visualizar los datos; en este caso se realiza a nivel estado, los cuales son 32 en México.

mx_state_map <- st_read("mexlatlong.shp")

## Reading layer `mexlatlong' from data source 
##   `C:\Users\AVRIL\Downloads\mexlatlong.shp' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -118.4042 ymin: 14.55055 xmax: -86.73862 ymax: 32.71846
## CRS:           NA

#mx_state_map <- read_sf("mexlatlong.shp")

Basic Map Making

Se realiza exploración visual básica del conjunto de datos geoespaciales a nivel estatal en México, útil para identificar patrones o validar la información unida mediante el geo_join. Cada mapa representa un atributo diferente, como:

• OBJECTID, FIPS_ADMIN, GMI_ADMIN, ADMIN_NAME: identificadores y nombres administrativos de los estados.

• POP_ADMIN: muestra diferencias de población entre estados (por ejemplo, CDMX y Edomex con colores más intensos).

• CNTRY_NAME, FIPS_CNTRY, GMI_CNTRY, TYPE_ENG: datos que son iguales para todos los estados, ya que todos pertenecen a México.

state_geodata <- geo_join(mx_state_map,mx_state,'OBJECTID','state_id',how='inner')
plot(state_geodata)

Visualización a nivel estado

Se observa y visualiza el mapa con los estados de México, acorde a esto se mezcla la data de las geocercas con alguna de las variables seleccionadas de la tabla de datos inicial para su análisis geoespacial.

tm_shape(mx_state_map) +
  tm_polygons(col = "black") +
  tm_compass(position = c("left", "bottom")) +
  tm_text("ADMIN_NAME", size = "AREA") +
  tm_layout(title = "Estados de México")

Distribución de Business Activity por Estado

businessactivity <- tm_shape(state_geodata) + 
  tm_polygons(
    col = "business_activity",
    palette = "-RdBu",  # el signo "-" invierte los colores para que los negativos estén en azul
    style = "cont",
    title = "Business Activity"
  ) +
  tm_layout(
    main.title = "Business Activity",
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

Distribución de Desempleo por Estado

unemployment <- tm_shape(state_geodata) + 
  tm_polygons(
    col = "unemployment",
    palette = "viridis",       # sin el "-" para que los valores más bajos sean claros
    style = "quantile",
    n = 8,
    title = "Unemployment"
  ) +
  tm_layout(
    main.title = "Unemployment",
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

Mapas de Business Activity & Unemployment

Business Activity:

• Casi todo el país está en azul, lo que sugiere que la mayoría de los estados tienen valores negativos o bajos de business_activity.

• Baja California aparece con un valor rojo oscuro → por lo que puede destacar en tener una actividad empresarial notablemente más alta o diferente.

Unemployment:

• Guerrero, Baja California, Campeche, tienen desempleo más bajo (morado oscuro).

• Nuevo León un rago medio-bajo de desempleo.

• Algunos estados como CDMX, Chiapas o Veracruz están en tonos más claros, lo que indica mayor desempleo relativo.

tmap_arrange(businessactivity, unemployment, ncol = 2)

Matriz de Conectividad Espacial

Se crea matriz de conectividad espacial o lista de vecinos entre polígonos para lo cual se usa el parámetro de QUEEN para definir la vecindad, en dónde se tomarán en cuaenta si dos polígonos son vecinos o, si comparten un borde o un punto (hace la cercanía más inclusiva).

Con base en el summary obtenemos lo siguiente:

• Hay 32 polígonos (32 estados de México).

• Hay 138 conexiones de vecindad entre ellos.

• Proporción de pares vecinos con relación al total posible = 13.5%

• Cada estado tiene, en promedio, 4.31 vecinos

• Sólo 1 estado tiene 1 vecino. Baja California Sur

• Sólo 1 estado está conectado con 9 vecinos San Luis Potosí

swm  <- poly2nb(mx_state_map, queen=T)

summary(swm) # The average number of neighbors is 4.31

## 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:
## 31 with 1 link
## 1 most connected region:
## 8 with 9 links

mx_state_map$ADMIN_NAME[8]   # Estado más conectado (con 9 vecinos)

## [1] "San Luis Potosi"

mx_state_map$ADMIN_NAME[31]  # Estado menos conectado (con 1 vecino)

## [1] "Baja California Sur"

state_geodata$ADMIN_NAME[8]   # Estado más conectado (con 9 vecinos)

## [1] "San Luis Potosi"

state_geodata$ADMIN_NAME[31]  # Estado menos conectado (con 1 vecino)

## [1] "Baja California Sur"

# Se muestran las 9 conexiones del Estado con mayor número de conexiones, es decir, de San Luis Potosí
swm[[8]]

## [1]  2  3  4  7  9 11 12 13 16

# Listado de estados conectados con San Luis Potosí
mx_state_map$ADMIN_NAME[swm[[8]]]

## [1] "Coahuila"   "Nuevo Leon" "Tamaulipas" "Zacatecas"  "Jalisco"   
## [6] "Veracruz"   "Guanajuato" "Queretaro"  "Hidalgo"

# Comprobación de estados conectados con San Luis Potosí
state_geodata$ADMIN_NAME[swm[[8]]]

## [1] "Coahuila"   "Nuevo Leon" "Tamaulipas" "Zacatecas"  "Jalisco"   
## [6] "Veracruz"   "Guanajuato" "Queretaro"  "Hidalgo"

Visualización Gráfica de la matriz de conectividad (Vecindades estilo ‘Reina’) del país.

# 1. Crear lista de vecinos
swm <- poly2nb(mx_state_map, queen = TRUE)

# 2. Convertir a matriz de pesos espaciales
sswm <- nb2listw(swm, style = "W", zero.policy = TRUE)

# 3. Convertir a objeto "Spatial" (para funciones más antigu
mx_state_map_a <- as(mx_state_map, "Spatial")

# 4. Obtener centroides (coordenadas)
mx_state_map_centroid <- coordinates(mx_state_map_a)

# 5. Graficar el mapa base
plot(mx_state_map_a, border = "blue", axes = FALSE, las = 1, main = "Mexico's States Queen SWM")
plot(mx_state_map_a, col = "grey", border = grey(0.9), axes = TRUE, add = TRUE)

# 6. Dibujar líneas de conectividad (vecindad)
plot(sswm, coords = mx_state_map_centroid, pch = 19, cex = 0.1, col = "red", add = TRUE)

# 7. Agregar nombres de los estados con ajuste
pointLabel(
  x = mx_state_map_centroid[,1],
  y = mx_state_map_centroid[,2],
  labels = mx_state_map$ADMIN_NAME,
  cex = 0.5,
  col = "black"
)

Análisis de Autocorrelación Espacial

Índice Global de Moran

Teóricamente el rango de Moran se muestra entre –1 y +1, en dónde: * +1 Alta autocorrelación espacial positiva (valores similares están cerca entre sí).

• 0 Distribución aleatoria (no hay patrón espacial).

• -1 Alta autocorrelación espacial negativa (valores diferentes están cerca entre sí).

Acorde a lo anterior, se genera codigo para determinar si el desempleo (unemployment) tiene autocorrelación espacial, es decir, si los estados ya sea con alto o bajo desempleo tienden a estar cerca unos de otros. También se analiza si el nivel de actividad empresarial (business_activity) tiene o no autocorrelación espacial, es decir, si los estados ya sea con alto o bajo nivel de actividad tienden a estar cerca unos de otros.

moran.test(state_geodata$business_activity, sswm)       # Global Moran's I is 0.18 but statistically significant (p-value > 10%). 

## 
##  Moran I test under randomisation
## 
## data:  state_geodata$business_activity  
## weights: sswm    
## 
## Moran I statistic standard deviate = 2.6024, p-value = 0.004629
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.184181268      -0.032258065       0.006917295

moran.test(state_geodata$unemployment, sswm)            # Global Moran's I is 0.021 but not statistically significant (p-value > 10%). 

## 
##  Moran I test under randomisation
## 
## data:  state_geodata$unemployment  
## weights: sswm    
## 
## Moran I statistic standard deviate = 0.43866, p-value = 0.3305
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.02162698       -0.03225806        0.01508982

Presentación de resultados

Actividad Empresarial (Business Activity )

• Moran’s I = 0.184 → hay autocorrelación espacial positiva: regiones con alta actividad empresarial tienden a estar cerca de otras similares.

• P-value = 0.0046 → es altamente significativo (p < 0.01).

Existe un patrón espacial fuerte. La actividad empresarial no está distribuida aleatoriamente; hay agrupamientos claros.

Desempleo (Unemployment)

• Moran’s I = 0.0216 (cercano a 0) sugiere poca autocorrelación espacial, regiones con alto desempleo pueden o no tender a estar cerca de otras similares, pero el índice no es igual o mayor a 0.35.

• P-value = 0.3305 > 0.10 → No es estadísticamente significativo (p>0.10).

No hay evidencia de que el desempleo esté agrupado geográficamente. Su distribución parece aleatoria pero no cumple con los estándares mínimos por ende, se realiza lag.

table <- data.frame(Variable = c("Business Activity", "Unemployment"), GM = c(0.184, 0.021), Significance = c("***", "NS"))
table

##            Variable    GM Significance
## 1 Business Activity 0.184          ***
## 2      Unemployment 0.021           NS

Identificación de Clusters

Se realiza cálculo de lags de acuerdo con la matriz de conectividad anteriormente creada

state_geodata$sp_lag_unemployment <- lag.listw(sswm, state_geodata$unemployment, zero.policy=TRUE) 
state_geodata$sp_lag_business_activity <- lag.listw(sswm, state_geodata$business_activity, zero.policy=TRUE)

Mapas de lag para Business Activity y Unemployment

Actividad Empresarial (Business Activity)

Mapa de businessactivity: Valor real por estado

• Baja California tiene actividad negativa más alta (rojo oscuro).

• El resto del país tiene valores más moderadamente negativos (azules claros a medios).

• Algunos estados del sur/sureste tienen actividad levemente positiva.

Mapa de ba_lag: Promedio de actividad empresarial en los estados vecinos

• El norte del país: colores fuertes en naranjas → se concentran los vecindarios (estados) con mayor actividad empresarial.

• Contraste con el sur: colores claros → sus vecinos tienen baja actividad empresarial.

• Algunos estados como Sonora y Baja California tienen valores altos en lag pero bajos en el presente, lo que sugiere que están rodeados de estados con mayor actividad.

Se observa que claramente si hay agrupamiento espacial positivo ya que se nota la diferencia de actividad empresarial entre el Norte y el Sur del país (clústeres regionales de actividad).

Desempleo (Unemployment)

Mapa de Unemployment: Porcentaje de desempleo en cada estado

• Distribución variada, no tan clara como Business Activity.

• Estados como Nayarit y Baja California Sur tienen desempleo relativamente alto.

• Varios estados del norte y centro tienen desempleo medio-bajo.

Mapa de Lag-1 Unemployment: Promedio de desempleo en los vecinos de cada estado

• Algunos estados como Chiapas o Campeche tienen valores moderados pero están rodeados de estados con bajo desempleo.

• Baja California Sur se mantiene con alto desempleo y sus vecinos también, reforzando la presencia de un cluster negativo.

• Muestra que el desempleo está más fragmentado, no hay patrones tan fuertes ni continuidad regional clara.

• Refuerza el resultado del estadístico Moran, que nos mostró que el desempleo no tiene autocorrelación espacial significativa, sino aleatoria.

Algunos estados están dentro de clusters donde tanto ellos como sus vecinos tienen valores similares, especialmente en desempleo (posible autocorrelación espacial). Asimismo, se observa que claramente no hay agrupamiento espacial positivo ya que no se diferencia la variación de la tasa de desempleo entre el Norte, Centro y Sur del país.

# Calcular rangos comunes
ba_range <- range(c(state_geodata$business_activity, state_geodata$sp_lag_business_activity), na.rm = TRUE)
unemp_range <- range(c(state_geodata$unemployment, state_geodata$sp_lag_unemployment), na.rm = TRUE)

# Mapa de Business Activity
businessactivity1 <- tm_shape(state_geodata) +  
  tm_polygons(
    col = "business_activity",
    palette = "-RdBu",
    style = "cont",
    title = "Business Activity",
    breaks = seq(ba_range[1], ba_range[2], length.out = 10)  # escala común
  ) +
  tm_layout(
    main.title = "Business Activity",
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mapa de Business Activity (Lag)
ba_lag1 <- tm_shape(state_geodata) + 
  tm_polygons(
    col = "sp_lag_business_activity",
    palette = "-RdBu",  # misma paleta
    style = "cont",
    title = "Business Activity (Lag-1)",
    breaks = seq(ba_range[1], ba_range[2], length.out = 10)  # misma escala
  ) +
  tm_layout(
    main.title = "Clusters of Business Activity (Lag-1)",
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mapa de Unemployment
unemployment1 <- tm_shape(state_geodata) + 
  tm_polygons(
    col = "unemployment",
    palette = "viridis",
    style = "cont",
    title = "Unemployment",
    breaks = seq(unemp_range[1], unemp_range[2], length.out = 8)  # escala común
  ) +
  tm_layout(
    main.title = "Unemployment",
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mapa de Unemployment (Lag)
unemployment_lag1 <- tm_shape(state_geodata) +  
  tm_polygons(
    col = "sp_lag_unemployment",
    palette = "viridis",  # misma paleta
    style = "cont",
    title = "Unemployment (Lag-1)",
    breaks = seq(unemp_range[1], unemp_range[2], length.out = 8)  # misma escala
  ) +
  tm_layout(
    main.title = "Clusters of Unemployment (Lag-1)",
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mostrar todos los mapas
tmap_arrange(businessactivity1, ba_lag1, unemployment1, unemployment_lag1, ncol = 2)

Scatterplots

Business Activity

Se realiza un modelo de regresión lineal, para identificar qué tan bien se puede predecir el valor de los estados vecinos con el valor del estado:

• Puntos: cada estado

• Eje X: actividad empresarial del estado

• Eje Y: el promedio de actividad empresarial de sus vecinos

• M1 (lm): La línea azul tiene pendiente positiva → hay agrupamiento espacial positivo, por lo que también refuerza el resultado de Moran que se arrojó previamente.

Con base al modelo, tiene una pendiente positiva, es decir, que los estados con más actividad tienden a estar cerca de otros similares → autocorrelación positiva.

M1 <- lm(sp_lag_business_activity ~ business_activity, state_geodata)

state_geodata$residuals <- abs(resid(M1))

n_outliers <- 5
top_outliers <- order(state_geodata$residuals, decreasing = TRUE)[1:n_outliers]

# Plot
plot(sp_lag_business_activity ~ business_activity, state_geodata, pch=21, asp=1, las=1, 
     col = "grey40", bg="grey80", main="Business Activity")

# Línea de regresión y líneas medias
abline(M1, col="blue")
abline(v = mean(state_geodata$business_activity), lty=3, col = "grey80")
abline(h = mean(state_geodata$sp_lag_business_activity), lty=3, col = "grey80")

# Puntos outliers
points(state_geodata$business_activity[top_outliers], 
       state_geodata$sp_lag_business_activity[top_outliers], 
       pch=21, bg="red", col="black", cex=1.4)

# Etiquetas con nombre (ajusta 'NAME' si tu columna tiene otro nombre)
text(state_geodata$business_activity[top_outliers],
     state_geodata$sp_lag_business_activity[top_outliers],
     labels = state_geodata$state[top_outliers],
     pos = 3, cex = 0.8, col = "black")

Unemployment

Se realiza un modelo de regresión lineal, para identificar como el desempleo en un estado predice el desempleo promedio de sus estados vecinos en donde:

• Puntos: cada estado

• Eje X: tasa de desempleo del estado

• Eje Y: el promedio de desempleo de sus vecinos

• M1 (lm): La línea azul tiene una pendiente casi plana (recta horizontalmente) → muy débil relación espacial, por lo que también refuerza el resultado de Moran que nos arrojó previamente

• Puntos dispersos → no hay agrupamientos claros. (No hay clústeres fuertes ni agrupamiento espacial de la variable.)

• El desempleo parece espacialmente aleatorio: no se agrupa de forma clara por regiones.

Con base al modelo, se determina que la variable tiene una pendiente muy baja → No hay evidencia de autocorrelación espacial significativa, es decir, un estado con alto desempleo no necesariamente está rodeado de otros con alto desempleo.

# Create a regression model
M2 <- lm(sp_lag_unemployment ~ unemployment, state_geodata)

# Calcular los residuos absolutos (distancia vertical a la línea de regresión)
state_geodata$residuals <- abs(resid(M2))

# Identificar los n puntos más lejanos (por ejemplo, los 5 mayores residuos)
n_outliers <- 5
top_outliers <- order(state_geodata$residuals, decreasing = TRUE)[1:n_outliers]

# Plot the data (Moran scatterplot) + LM
plot(sp_lag_unemployment ~ unemployment, state_geodata, pch=21, asp=1, las=1, 
     col = "grey40", bg="grey80", main="Unemployment")
abline(M1, col="blue") # Línea de regresión
abline(v = mean(state_geodata$unemployment), lty=3, col = "grey80")
abline(h = mean(state_geodata$sp_lag_unemployment), lty=3, col = "grey80")

points(state_geodata$unemployment[top_outliers], 
       state_geodata$sp_lag_unemployment[top_outliers], 
       pch=21, bg="red", col="black", cex=1.4)

# Etiquetas con nombre (ajusta 'NAME' si tu columna tiene otro nombre)
text(state_geodata$unemployment[top_outliers],
     state_geodata$sp_lag_unemployment[top_outliers],
     labels = state_geodata$state[top_outliers],
     pos = 3, cex = 0.8, col = "black")

Visualización Espacial de Clusters (Conglomerados)

Segregación en las diferentes posibilidades de agrupaciones:

• No significante

• HotSpots: High - High

• ColdSpots: Low - Low

• Atípicos: High - Low

• Atípicos: Low - High

Business Activity

Se generan clústers:

• Low-High: Baja actividad empresarial, pero rodeado de vecinos con altactividad empresarial (posible outlier espacial).

• High-High: Alta actividad empresarial, rodeado de vecinos con alta actividad empresarial (clúster positivo).

• High-Low: Alta actividad empresarial, pero vecinos con baja actividad empresarial (posible outlier espacial).

• Low-Low: Baja actividad empresarial, rodeado de vecinos con baja actividad empresarial (clúster negativo).

Con base a la colorimetría del mapa se determina que:

• 🔴 Estados sin evidencia de agrupamiento espacial

• 🟢 Sonora estado con alta actividad empresarial rodeados de otros similares

• 🔵 Península de Yucatán estados con baja actividad empresarial rodeado de otros similares (clúster negativo)

• 🟣 Baja California Sur estado con baja actividad empresarial rodeado de vecinos con alta actividad (outlier negativo)

swm_a <- queen_weights(mx_state_map) # queen spatial weight matrix (alternative format)

# 1. Calcular LISA
lisa_business <- local_moran(swm_a, state_geodata["business_activity"])

# 2. Asignar clústeres con etiquetas
state_geodata$cluster_business <- as.factor(lisa_business$GetClusterIndicators())
levels(state_geodata$cluster_business) <- lisa_business$GetLabels()

# 3. Mapa con estilo por defecto
library(ggplot2)

ggplot(data = state_geodata) +
  geom_sf(aes(fill = cluster_business), color = "white") +
  ggtitle(label = "Business Activity", subtitle = "Mexico's States") +
  theme_minimal()

Unemployment

Se generan clústers:

• Low-High: Bajo desempleo, pero rodeado de vecinos con alto desempleo (posible outlier espacial).

• High-High: Alto desempleo, rodeado de vecinos con alto desempleo (clúster positivo).

• High-Low: Alto desempleo, pero vecinos con bajo desempleo (posible outlier espacial).

• Low-Low: Bajo desempleo, rodeado de vecinos con bajo desempleo (clúster negativo).

Con base a la colorimetría del mapa se determina que:

• 🔴 Estados sin evidencia de agrupamiento espacial

• 🟢 Hidalgo estado con alto desempleo rodeados de otros similares

• 🔵 Quintana Roo estado con baja tasa de desempleo rodeado de otros similares (clúster negativo)

sswm_a <- queen_weights(mx_state_map) # queen spatial weight matrix (alternative format)
lisa_income <- local_moran(sswm_a, state_geodata["unemployment"]) 
state_geodata$cluster_unemployment <- as.factor(lisa_income$GetClusterIndicators())
levels(state_geodata$cluster_unemployment)<-lisa_income$GetLabels() 

ggplot(data=state_geodata) +
  geom_sf(aes(fill=cluster_unemployment)) + 
  ggtitle(label = "Unemployment", subtitle = "Mexico's States")

Modelos

Se cargan librerías para el ajuste de modelos de regresión:

library(spatialreg)
library(stargazer)

Regresión Simple

Se busca explicar la variable de Inversión Extranjera Directa considerando las variables de Actividad Empresarial y la Tasa de Desempleo.

El modelo_a, aunque significativo, realmente no es tan bueno, pues la \(R^2\) ajustada apenas sobrepasa el 0.5 (\(\bar{R^2}\) = 0.5821).

Con base a esto, se generaán otros modelos para dar una mayor complejidad al considerar terminos de autocorrelación espacial.

model_a <- lm(new_fdi_real_mxn ~ business_activity + unemployment, data = state_geodata)
summary(model_a)

## 
## Call:
## lm(formula = new_fdi_real_mxn ~ business_activity + unemployment, 
##     data = state_geodata)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -8473  -4523  -2443   1047  26134 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)          12176       4971   2.449   0.0206 *
## business_activity     4695       1722   2.727   0.0107 *
## unemployment         97524     124597   0.783   0.4401  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8077 on 29 degrees of freedom
## Multiple R-squared:  0.2043, Adjusted R-squared:  0.1494 
## F-statistic: 3.722 on 2 and 29 DF,  p-value: 0.03641

AIC(model_a)

## [1] 671.4531

SAR - Spatial AutoRegressive Model

Se realiza modelo SAR, también conocido como Spatial Lag Model Esencial, en el cual se incluye la matriz de conectividad espacial dada en el parámetro ‘listw’ por el objeto previamente calculado con el nombre de ‘sswm’. Se observa que el AIC es mayor que el de regresión simple.

model_b <- lagsarlm(new_fdi_real_mxn ~ business_activity + unemployment, data = state_geodata, listw = sswm) 
summary(model_b)

## 
## Call:lagsarlm(formula = new_fdi_real_mxn ~ business_activity + unemployment, 
##     data = state_geodata, listw = sswm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -8378.93 -4533.51 -2475.15   922.67 26266.71 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)        11811.5     5539.3  2.1323 0.032982
## business_activity   4623.7     1740.1  2.6572 0.007879
## unemployment       97513.3   118721.1  0.8214 0.411438
## 
## Rho: 0.032956, LR test value: 0.017546, p-value: 0.89462
## Approximate (numerical Hessian) standard error: 0.26295
##     z-value: 0.12534, p-value: 0.90026
## Wald statistic: 0.015709, p-value: 0.90026
## 
## Log likelihood: -331.7178 for lag model
## ML residual variance (sigma squared): 59069000, (sigma: 7685.6)
## Number of observations: 32 
## Number of parameters estimated: 5 
## AIC: 673.44, (AIC for lm: 671.45)

AIC(model_b)

## [1] 673.4356

SEM - Spatial Error Model

Se genera Modelo de Error Espacial. Igual que en el modelo SAR es esencial, incluir la matriz de conectividad espacial dada en el parámetro ‘listw’ por el objeto previamente calculado con el nombre de ‘sswm’. Se calcula que su valor de AIC es menor que el SAR y con una diferencia mayor de 0.21 al AIC de la regresión simple.

model_c <- errorsarlm(new_fdi_real_mxn ~ business_activity + unemployment, data = state_geodata, listw = sswm)
summary(model_c)

## 
## Call:errorsarlm(formula = new_fdi_real_mxn ~ business_activity + unemployment, 
##     data = state_geodata, listw = sswm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9664.0 -4067.2 -2083.9  1860.1 25175.2 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                   Estimate Std. Error z value  Pr(>|z|)
## (Intercept)        13904.6     4153.0  3.3481 0.0008137
## business_activity   5786.5     1386.1  4.1746 2.985e-05
## unemployment      100452.6   110678.7  0.9076 0.3640865
## 
## Lambda: -0.39762, LR test value: 1.7915, p-value: 0.18074
## Asymptotic standard error: 0.23953
##     z-value: -1.66, p-value: 0.096917
## Wald statistic: 2.7556, p-value: 0.096917
## 
## Log likelihood: -330.8308 for error model
## ML residual variance (sigma squared): 53778000, (sigma: 7333.3)
## Number of observations: 32 
## Number of parameters estimated: 5 
## AIC: 671.66, (AIC for lm: 671.45)

AIC(model_c)

## [1] 671.6616

SDM - Spatial Durbin Model

Igual que en los modelo SAR y SEM es esencial, incluir la matriz de conectividad espacial dada en el parámetro ‘listw’ por el objeto previamente calculado con el nombre de ‘sswm’ y aquí en particular se debe añadir el parámetro ‘type’ con valor ‘mixed’ para que en esencia mezcle los dos modelos previos (SAR y SEM). SU valor de AIC es el más alto pero, logra la menor varianza residual y mejor interpretación estructural.

model_d <- lagsarlm(new_fdi_real_mxn ~ business_activity + unemployment, data = state_geodata, listw = sswm, type="mixed") 
summary(model_d)

## 
## Call:lagsarlm(formula = new_fdi_real_mxn ~ business_activity + unemployment, 
##     data = state_geodata, listw = sswm, type = "mixed")
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9395.9 -4069.3 -2000.3  2234.7 25838.6 
## 
## Type: mixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                       Estimate Std. Error z value Pr(>|z|)
## (Intercept)            16582.9    10206.0  1.6248 0.104201
## business_activity       4471.6     1664.1  2.6871 0.007208
## unemployment          105464.9   110826.1  0.9516 0.341287
## lag.business_activity   4018.2     2206.8  1.8208 0.068638
## lag.unemployment       97809.4   261432.2  0.3741 0.708308
## 
## Rho: -0.23256, LR test value: 0.6703, p-value: 0.41295
## Approximate (numerical Hessian) standard error: 0.28212
##     z-value: -0.82434, p-value: 0.40974
## Wald statistic: 0.67954, p-value: 0.40974
## 
## Log likelihood: -329.7602 for mixed model
## ML residual variance (sigma squared): 51591000, (sigma: 7182.7)
## Number of observations: 32 
## Number of parameters estimated: 7 
## AIC: 673.52, (AIC for lm: 672.19)

AIC(model_d)

## [1] 673.5205

Comparativa

Aunque el modelo c presenta un mejor rendimiento estadístico con el menor AIC (671.66), el objetivo principal es identificar cómo el entorno —específicamente la actividad de los estados vecinos— influye en la Inversión Extranjera Directa (FDI). Por ello, se selecciona el modelo d, ya que permite captar efectos espaciales indirectos (influencia de regiones vecinas), a pesar de tener un AIC ligeramente superior. Además, tanto la variable business_activity como su retardo espacial resultan significativas o cercanas a la significancia, lo cual evidencia una interdependencia regional en la atracción de inversión. De esta manera, el modelo d ofrece una representación más completa y realista del fenómeno espacial analizado.

stargazer(model_a, model_b, model_c, model_d, type = "text", title="Estimated Regression Results")

## 
## Estimated Regression Results
## =======================================================================================
##                                              Dependent variable:                       
##                       -----------------------------------------------------------------
##                                               new_fdi_real_mxn                         
##                               OLS             spatial        spatial        spatial    
##                                            autoregressive     error      autoregressive
##                               (1)               (2)            (3)            (4)      
## ---------------------------------------------------------------------------------------
## business_activity         4,695.418**       4,623.750***   5,786.464***   4,471.587*** 
##                           (1,722.100)       (1,740.075)    (1,386.100)    (1,664.106)  
##                                                                                        
## unemployment               97,523.970        97,513.350    100,452.600    105,464.900  
##                          (124,597.400)     (118,721.100)  (110,678.700)  (110,826.100) 
##                                                                                        
## lag.business_activity                                                      4,018.200*  
##                                                                           (2,206.836)  
##                                                                                        
## lag.unemployment                                                           97,809.420  
##                                                                          (261,432.200) 
##                                                                                        
## Constant                  12,175.640**      11,811.500**  13,904.580***    16,582.930  
##                           (4,971.108)       (5,539.305)    (4,152.975)    (10,206.020) 
##                                                                                        
## ---------------------------------------------------------------------------------------
## Observations                   32                32             32             32      
## R2                           0.204                                                     
## Adjusted R2                  0.149                                                     
## Log Likelihood                                -331.718       -330.831       -329.760   
## sigma2                                     59,068,556.000 53,777,605.000 51,590,758.000
## Akaike Inf. Crit.                             673.436        671.662        673.520    
## Residual Std. Error   8,076.685 (df = 29)                                              
## F Statistic           3.722** (df = 2; 29)                                             
## Wald Test (df = 1)                             0.016          2.756*         0.680     
## LR Test (df = 1)                               0.018          1.792          0.670     
## =======================================================================================
## Note:                                                       *p<0.1; **p<0.05; ***p<0.01

Conclusión

Se evaluaron cuatro modelos de regresión (OLS y espaciales) para explicar la variación de la inversión extranjera directa (new_fdi_real_mxn) en los estados de México. El modelo d, que incorpora efectos espaciales autoregresivos y los lags espaciales de los predictores, fue seleccionado como el más adecuado por las siguientes razones:

• Captura la dependencia espacial de la FDI: incluye un componente autoregresivo que reconoce que la inversión en un estado puede estar influida por la de sus estados vecinos.

• Incorpora los efectos del entorno: al incluir el lag.business_activity, considera que la actividad empresarial en los estados vecinos también afecta las decisiones de inversión.

• Mantenimiento del ajuste estadístico: aunque el modelo c tuvo el AIC más bajo, el modelo d conserva un ajuste competitivo (AIC = 673), con todos los efectos relevantes incluidos.

• Interpretación más amplia y adecuada al entoro de estudio: permite entender no solo los efectos directos, sino también los efectos espaciales indirectos, lo que es clave para decisiones de política económica regional.

Referencias

What is exploratory data analysis? https://www.ibm.com/topics/exploratory-data-analysis

Exploratory Spatial Data Analysis https://cran.r-project.org/web/packages/geostan/vignettes/measuring-sa.html