library(foreign)                  # import external files
library(dplyr)                    # data manipulation 
library(spdep)                    # a collection of functions to create spatial weight matrix 
library(tigris)                   # allows to work with shapeflies
library(rgeoda)                   # spatial data analysis based on GeoDa
library(RColorBrewer)             # offers several color palettes 
library(viridis)                  # offers several color palettes
library(ggplot2)                  # to create plots and graphics from dataset 
library(tmap)                     # making maps so spatial data distributions are visualized
library(sf)                       # functions to encode spatial vector data 
library(sp)                       # classes and methods for spatial data
library(spgwr)                    # geographically weighted regresion
library(GWmodel)                  # require for kernel adaptive creation to apply GWR
mx_state <- read.csv("state_data.csv")
mx_mpio  <- read.csv("mpio_data.csv")
mx_state_map <- st_read("mexlatlong.shp")
## Reading layer `mexlatlong' from data source 
##   `C:\Users\Alberto JIménez\Downloads\Modulo Rodolfo\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 84
mx_mpio_map <- st_read("Mexican Municipalities.shp")
## Reading layer `Mexican Municipalities' from data source 
##   `C:\Users\Alberto JIménez\Downloads\Modulo Rodolfo\Mexican Municipalities.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2456 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -118.4076 ymin: 14.5321 xmax: -86.71041 ymax: 32.71865
## Geodetic CRS:  WGS 84
state_geodata <- geo_join(mx_state_map,mx_state,'OBJECTID','state_id',how='inner')
mpio_geodata  <- geo_join(mx_mpio_map,mx_mpio,'IDUNICO','clave_municipio',how='inner')

Distribución de Densidad Empleo por Estado

layout_clean <- tm_layout(
  legend.outside = TRUE,
  legend.outside.position = "left",
  title.position = c('right', 'top'),
  title.size = 1
)

pop_emp <- tm_shape(state_geodata) + 
  tm_polygons(col = "employment", palette="Blues", style="quantile", 
              n=8, title="Employment Density") +
  layout_clean +
  tm_layout(main.title = 'Employment Density')

Distribución de Densidad Salarial por Estado

pop_wage <- tm_shape(state_geodata) + 
  tm_polygons(col = "real_wage", palette="-viridis", style="quantile", 
              n=8, title="Wage Density") +
  layout_clean +
  tm_layout(main.title = 'Wage Density')
tmap_arrange(pop_emp, pop_wage, ncol = 2)

Matriz de Conectividad Espacial

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

Visualización Gráfica de la matriz de conectividad (Vecindades estilo ‘Reina’)

sswm <- nb2listw(swm, style="W", zero.policy = TRUE)

mx_state_map_a <- as(mx_state_map, "Spatial")
mx_state_map_centroid <- coordinates(mx_state_map_a) 
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=T,add=T) 
plot(sswm,coords=mx_state_map_centroid,pch=19,cex=0.1,col="red",add=T)

Índice Globarl de Moran

moran.test(state_geodata$employment, sswm)            
## 
##  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.01474416
moran.test(state_geodata$real_wage, sswm)                 
## 
##  Moran I test under randomisation
## 
## data:  state_geodata$real_wage  
## weights: sswm    
## 
## Moran I statistic standard deviate = 0.86376, p-value = 0.1939
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.07246688       -0.03225806        0.01469995
table <- data.frame(Variable = c("Employment", "Real_wage"),GM = c(0.0202, 0.0724),Significance = c("NS", "NS"))
table
##     Variable     GM Significance
## 1 Employment 0.0202           NS
## 2  Real_wage 0.0724           NS
state_geodata$Lag_Employment <- lag.listw(sswm, state_geodata$employment, zero.policy=TRUE) 
state_geodata$Lag_Real_Wage  <- lag.listw(sswm, state_geodata$real_wage, zero.policy=TRUE)

Mapas de lag para el Empleo y Salario

layout_clean <- tm_layout(legend.outside = TRUE,
                           legend.outside.position = "left",
                           title.position = c('right', 'top'),
                           title.size = 1)

pop_emp_lag <- tm_shape(state_geodata) + 
  tm_polygons(col = "Lag_Employment", palette="Blues", 
              style="quantile", n=8, title="Employment") + layout_clean +
  tm_layout(main.title = 'Clusters of Employment')

pop_wage_lag <- tm_shape(state_geodata) + 
  tm_polygons(col = "Lag_Real_Wage", palette="-viridis", 
              style="quantile", n=8, title="Real Wage") + layout_clean +
  tm_layout(main.title = 'Clusters of Real Wage')

tmap_arrange(pop_emp, pop_emp_lag, pop_wage, pop_wage_lag, ncol = 2)

Employment

M_emp <- lm(Lag_Employment ~ employment, data = state_geodata)

# Gráfico de Dispersión
plot(Lag_Employment ~ employment, data = state_geodata, 
     pch=21, asp=1, las=1, col = "grey40", bg="grey80", 
     main="Moran Scatterplot: Empleo",
     xlab="Empleo (Original)", ylab="Lag de Empleo")

# Líneas de referencia
abline(M_emp, col="blue", lwd=2) # Línea de regresión
abline(v = mean(state_geodata$employment), lty=3, col = "red") # Media X
abline(h = mean(state_geodata$Lag_Employment), lty=3, col = "red") # Media Y

M_wage <- lm(Lag_Real_Wage ~ real_wage, data = state_geodata)

# Gráfico de Dispersión
plot(Lag_Real_Wage ~ real_wage, data = state_geodata, 
     pch=21, asp=1, las=1, col = "grey40", bg="grey80", 
     main="Moran Scatterplot: Salario Real",
     xlab="Salario Real (Original)", ylab="Lag de Salario")

# Líneas de referencia
abline(M_wage, col="darkgreen", lwd=2) # Línea de regresión
abline(v = mean(state_geodata$real_wage), lty=3, col = "red") # Media X
abline(h = mean(state_geodata$Lag_Real_Wage), lty=3, col = "red") # Media Y

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

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

sswm_a <- queen_weights(mx_state_map) 
lisa_wage <- local_moran(sswm_a, state_geodata["real_wage"]) 
state_geodata$cluster_wage <- as.factor(lisa_wage$GetClusterIndicators())
levels(state_geodata$cluster_wage)<-lisa_wage$GetLabels() 

ggplot(data=state_geodata) +
  geom_sf(aes(fill=cluster_wage)) + 
  ggtitle(label = "Average Wage ", subtitle = "Mexico's States")