Actividad 02 - Análisis Espacial
Base de datos: state_data.csv + mexlatlong.shp… (variable a predecir: Employment y Log_new_fdi_real_mxn)
EQUIPO 1
2026-10-06
1. Objetivo
Analizar la distribución espacial y la autocorrelación espacial de dos variables a nivel estado en México:
- employment: tasa de empleo.
- log_new_fdi_real_mxn: logaritmo de los nuevos flujos de inversión extranjera directa (IED) en pesos reales.
2. Librerías
Este bloque instala automáticamente solo los paquetes que falten (con el repositorio CRAN ya indicado) y luego los carga:
paquetes <- c("dplyr", "sf", "spdep", "tigris", "ggplot2", "viridis", "gridExtra")
faltan <- paquetes[!paquetes %in% rownames(installed.packages())]
if (length(faltan) > 0) install.packages(faltan, repos = "https://cloud.r-project.org")
library(dplyr)
library(sf)
library(spdep)
library(tigris) # geo_join
library(ggplot2)
library(viridis)
library(gridExtra)3. Datos
Los archivos state_data.csv, mexlatlong.shp
(y sus archivos acompañantes .dbf, .shx, .prj, .cpg) deben estar en la
misma carpeta que este .Rmd.
# Carpeta donde están el CSV y los archivos mexlatlong.*
ruta <- "C:/TEC 8 semestre/"
# Verificación: debe listar al menos .shp, .shx y .dbf
list.files(ruta, pattern = "mexlatlong")## [1] "mexlatlong.cpg" "mexlatlong.dbf" "mexlatlong.prj"
## [4] "mexlatlong.sbn" "mexlatlong.sbx" "mexlatlong.shp"
## [7] "mexlatlong.shx" "mexlatlong_filesss" "mexlatlong_filesss.zip"mx_state <- read.csv(paste0(ruta, "state_data(in).csv"))
mx_state_map <- st_read(paste0(ruta, "mexlatlong.shp"))## Reading layer `mexlatlong' from data source `C:\TEC 8 semestre\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
## Geodetic CRS: WGS 84glimpse(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, …glimpse(mx_state_map)## Rows: 32
## Columns: 20
## $ OBJECTID <dbl> 888, 933, 976, 978, 998, 1004, 1026, 1034, 1051, 1057, 1058…
## $ FIPS_ADMIN <chr> "MX06", "MX07", "MX19", "MX28", "MX25", "MX10", "MX32", "MX…
## $ GMI_ADMIN <chr> "MEX-CHH", "MEX-CDZ", "MEX-NLE", "MEX-TML", "MEX-SIN", "MEX…
## $ ADMIN_NAME <chr> "Chihuahua", "Coahuila", "Nuevo Leon", "Tamaulipas", "Sinal…
## $ FIPS_CNTRY <chr> "MX", "MX", "MX", "MX", "MX", "MX", "MX", "MX", "MX", "MX",…
## $ GMI_CNTRY <chr> "MEX", "MEX", "MEX", "MEX", "MEX", "MEX", "MEX", "MEX", "ME…
## $ CNTRY_NAME <chr> "Mexico", "Mexico", "Mexico", "Mexico", "Mexico", "Mexico",…
## $ POP_ADMIN <dbl> 2656214, 2145539, 3370912, 2272724, 2397706, 1467826, 13836…
## $ TYPE_ENG <chr> "State", "State", "State", "State", "State", "State", "Stat…
## $ TYPE_LOC <chr> "Estado", "Estado", "Estado", "Estado", "Estado", "Estado",…
## $ SQKM <dbl> 247935.02, 150843.95, 65173.05, 79502.24, 57638.85, 120674.…
## $ SQMI <dbl> 95727.70, 58240.87, 25163.31, 30695.81, 22254.36, 46592.46,…
## $ COLOR_MAP <chr> "12", "2", "3", "11", "5", "4", "7", "1", "8", "4", "9", "1…
## $ Shape_Leng <dbl> 22.609277, 18.993090, 15.426171, 18.023144, 16.466051, 17.5…
## $ Shape_Area <dbl> 22.8909849, 13.7336548, 5.8446681, 7.0565626, 5.1455245, 10…
## $ OBJECTID_1 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
## $ OBJECTID_2 <dbl> 888, 933, 976, 978, 998, 1004, 1026, 1034, 1051, 1057, 1058…
## $ longitude <dbl> -106.44431, -102.03356, -99.83125, -98.62181, -107.48280, -…
## $ latitude <dbl> 28.80303, 27.30662, 25.60105, 24.29819, 25.02179, 24.93519,…
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-103.6309 2..., MULTIPOLYGON (…Unión de la tabla con las geocercas (misma llave que en clase: OBJECTID y state_id):
state_geodata <- geo_join(mx_state_map, mx_state, "OBJECTID", "state_id", how = "inner")
nrow(state_geodata) # debe ser 32## [1] 32sum(is.na(state_geodata$employment)) # debe ser 0## [1] 0Resumen de las variables
summary(state_geodata[, c("employment", "log_new_fdi_real_mxn")] %>% st_drop_geometry())## employment log_new_fdi_real_mxn
## Min. :0.9400 Min. :-2.680
## 1st Qu.:0.9600 1st Qu.: 3.152
## Median :0.9700 Median : 3.665
## Mean :0.9684 Mean : 3.320
## 3rd Qu.:0.9725 3rd Qu.: 3.980
## Max. :0.9900 Max. : 4.5404. Matriz de conectividad espacial (vecindad tipo Reina)
swm <- poly2nb(state_geodata, queen = TRUE)
summary(swm)## 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 linkssswm <- nb2listw(swm, style = "W", zero.policy = TRUE)centroides <- st_coordinates(st_centroid(st_geometry(state_geodata)))
plot(st_geometry(state_geodata), col = "grey90", border = "white",
main = "Matriz de vecindad Reina - Estados de México")
plot(swm, coords = centroides, add = TRUE, col = "red", pch = 19, cex = 0.4)
5. Rezago espacial (lag-1)
El lag-1 de cada estado es el promedio ponderado del valor de sus vecinos.
state_geodata$lag_employment <- lag.listw(sswm, state_geodata$employment, zero.policy = TRUE)
state_geodata$lag_log_fdi <- lag.listw(sswm, state_geodata$log_new_fdi_real_mxn, zero.policy = TRUE)6. Mapas
mapa <- function(var, titulo, leyenda, opcion) {
ggplot(state_geodata) +
geom_sf(aes(fill = .data[[var]]), color = "white", linewidth = 0.2) +
scale_fill_viridis_c(option = opcion, direction = -1, name = leyenda) +
labs(title = titulo) +
theme_minimal() +
theme(legend.position = "bottom")
}
m1 <- mapa("employment", "Employment (original)", "Employment", "mako")
m2 <- mapa("lag_employment", "Employment (lag-1 espacial)", "Lag Employment", "mako")
grid.arrange(m1, m2, ncol = 2)
m3 <- mapa("log_new_fdi_real_mxn", "Log Nueva IED real (original)", "Log IED", "viridis")
m4 <- mapa("lag_log_fdi", "Log Nueva IED real (lag-1 espacial)", "Lag Log IED", "viridis")
grid.arrange(m3, m4, ncol = 2)
7. Prueba global de Moran
H0: no hay autocorrelación espacial (distribución espacial aleatoria). H1: sí hay autocorrelación espacial.
Employment
moran_emp <- moran.test(state_geodata$employment, sswm, zero.policy = TRUE)
moran_emp##
## Moran I test under randomisation
##
## data: state_geodata$employment
## weights: sswm
##
## Moran I statistic standard deviate = 0.43241, p-value = 0.3327
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.02024763 -0.03225806 0.01474416Log_new_fdi_real_mxn
moran_fdi <- moran.test(state_geodata$log_new_fdi_real_mxn, sswm, zero.policy = TRUE)
moran_fdi##
## Moran I test under randomisation
##
## data: state_geodata$log_new_fdi_real_mxn
## weights: sswm
##
## Moran I statistic standard deviate = 0.15067, p-value = 0.4401
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.018477521 -0.032258065 0.008365411Tabla resumen
resumen <- data.frame(
Variable = c("Employment", "Log_new_fdi_real_mxn"),
Moran_I = c(moran_emp$estimate[["Moran I statistic"]],
moran_fdi$estimate[["Moran I statistic"]]),
Esperado = c(moran_emp$estimate[["Expectation"]],
moran_fdi$estimate[["Expectation"]]),
p_value = c(moran_emp$p.value, moran_fdi$p.value)
)
resumen$Significancia <- cut(resumen$p_value,
breaks = c(-Inf, 0.01, 0.05, 0.10, Inf),
labels = c("***", "**", "*", "NS"))
resumenPrueba de permutaciones (Monte Carlo), como verificación de robustez
set.seed(123)
moran.mc(state_geodata$employment, sswm, nsim = 999, zero.policy = TRUE)##
## Monte-Carlo simulation of Moran I
##
## data: state_geodata$employment
## weights: sswm
## number of simulations + 1: 1000
##
## statistic = 0.020248, observed rank = 696, p-value = 0.304
## alternative hypothesis: greatermoran.mc(state_geodata$log_new_fdi_real_mxn, sswm, nsim = 999, zero.policy = TRUE)##
## Monte-Carlo simulation of Moran I
##
## data: state_geodata$log_new_fdi_real_mxn
## weights: sswm
## number of simulations + 1: 1000
##
## statistic = -0.018478, observed rank = 571, p-value = 0.429
## alternative hypothesis: greater8. Diagramas de dispersión de Moran
moran.plot(state_geodata$employment, listw = sswm, zero.policy = TRUE,
labels = state_geodata$state, pch = 19,
xlab = "Employment", ylab = "Lag-1 de Employment",
main = "Moran Scatterplot - Employment")
moran.plot(state_geodata$log_new_fdi_real_mxn, listw = sswm, zero.policy = TRUE,
labels = state_geodata$state, pch = 19,
xlab = "Log Nueva IED real", ylab = "Lag-1 de Log Nueva IED real",
main = "Moran Scatterplot - Log IED")
Versión alternativa en ggplot2 (valores estandarizados, con cuadrantes):
scatter <- function(x, lagx, titulo) {
d <- data.frame(z = as.numeric(scale(x)), lagz = as.numeric(scale(lagx)),
estado = state_geodata$state)
ggplot(d, aes(z, lagz, label = estado)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_point(color = "steelblue", size = 2) +
geom_smooth(method = "lm", se = FALSE, color = "red", linewidth = 0.7) +
geom_text(size = 2.5, vjust = -0.8, check_overlap = TRUE) +
labs(title = titulo, x = "Valor estandarizado", y = "Lag-1 estandarizado") +
theme_minimal()
}
grid.arrange(
scatter(state_geodata$employment, state_geodata$lag_employment, "Employment"),
scatter(state_geodata$log_new_fdi_real_mxn, state_geodata$lag_log_fdi, "Log Nueva IED real"),
ncol = 2)
9. Interpretación (completar después de ejecutar)
- Employment: I de Moran = , p-valor = -> ___ (autocorrelación positiva / negativa / no significativa). Estados o regiones que forman clusters: ___
- Log_new_fdi_real_mxn: I de Moran = , p-valor = -> . Estados o regiones que forman clusters:
- Nota sobre outliers: Durango tiene un valor atípico muy bajo de log IED (-2.68), que puede afectar el estadístico de Moran y el scatterplot; conviene mencionarlo y, si se desea, repetir el análisis sin ese estado como prueba de sensibilidad.