[Descargar - RMarkdown]

Exploración de datos

Leer los paquetes estadísticos:

library(sp)
library(sf) 
library(spData)
library(spdep)
library(spatialreg)

#Gráficas
library(ggplot2)
library(tmap) # Mapas estaticos e interactivos
library(mapview) # bueno para explorar pero dificil ajustar

Primero usar la libreria de sf para mapear y explorar nuestros datos:

muni = st_read("./data/mx_city_metro.shp")
## Reading layer `mx_city_metro' from data source 
##   `/Users/irenefarah/Downloads/sp_econ/data/mx_city_metro.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 76 features and 40 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2745632 ymin: 774927.1 xmax: 2855437 ymax: 899488.5
## proj4string:   +proj=lcc +lat_1=17.5 +lat_2=29.5 +lat_0=12 +lon_0=-102 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs

Pueden usar attach(muni) para referirse a las variables por nombre, pero no lo haré en este caso.

Formato puede ser punto, línea, polígono.

Datos provienen de:

Ver proyecciones de los datos:

st_crs(muni)
## Coordinate Reference System:
##   No user input
##   wkt:
## PROJCS["Lambert_Conformal_Conic",
##     GEOGCS["GCS_GRS 1980(IUGG, 1980)",
##         DATUM["unknown",
##             SPHEROID["GRS80",6378137,298.257222101]],
##         PRIMEM["Greenwich",0],
##         UNIT["Degree",0.017453292519943295]],
##     PROJECTION["Lambert_Conformal_Conic_2SP"],
##     PARAMETER["standard_parallel_1",17.5],
##     PARAMETER["standard_parallel_2",29.5],
##     PARAMETER["latitude_of_origin",12],
##     PARAMETER["central_meridian",-102],
##     PARAMETER["false_easting",2500000],
##     PARAMETER["false_northing",0],
##     UNIT["Meter",1]]

Transformación de la proyección:

municipios_4326 = st_transform(muni, 4326)
st_crs(municipios_4326)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

Crear mapas coropléticos de distintas formas:

ggplot(municipios_4326, aes(fill = ppov_20)) + 
  geom_sf() +  # tells ggplot that geographic data are being plotted
  scale_fill_viridis_c() +
  theme_minimal() + 
  labs(title = "Porcentaje de pobreza, 2020")

# Mapa interactivo
mapview(municipios_4326, zcol = 'ppov_20')
# Agregar los dos
mapview(municipios_4326, zcol = 'ppov_20') +
  mapview(muni, color = 'red', alpha.regions = 0, legend=FALSE)
#tmap_mode('plot') # estaticos
tmap_mode('view') # interactivo
## tmap mode set to interactive viewing
tm_shape(municipios_4326) + 
  tm_polygons(col = 'tan', 
              border.col = 'darkgreen', 
              alpha = 0.5) # transparencia

Guía de código de colores en R:

Examinar nuestros datos:

tmap_mode('plot') # estaticos

pe = tm_shape(municipios_4326) + 
  tm_polygons(col = 'pepov_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Pobreza extrema, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

ping1 = tm_shape(municipios_4326) + 
  tm_polygons(col = 'ping1_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Por debajo línea ingresos, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

educ = tm_shape(municipios_4326) + 
  tm_polygons(col = 'peduc_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Educación, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)
  
salud = tm_shape(municipios_4326) + 
  tm_polygons(col = 'pheal_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Salud, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

sbv = tm_shape(municipios_4326) + 
  tm_polygons(col = 'pqdwel_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Servicios básicos de la vivienda, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

alim = tm_shape(municipios_4326) + 
  tm_polygons(col = 'pfood_20', 
              border.col = 'darkgreen', 
              alpha = 0.5) + # transparencia
    tm_layout("Inseguridad Alimentaria, 2020",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

tmap_arrange(pe, ping1, educ, salud, sbv, alim, ncol=3, nrow=2)

Pesos espaciales en R

Leer pesos espaciales de archivos existentes (.gal y .gwt):

#Importar de archivo existente (abrir archivo)
#Pasar de "neighbor lists" a "weight lists" (listw object)
queen.nb_existente = read.gal("./data/mx_city_metro_q1.gal", region.id=muni$ID,override.id=TRUE) #spdep

summary(queen.nb_existente)
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 400 
## Percentage nonzero weights: 6.925208 
## Average number of links: 5.263158 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  2  3 11 13 16 12  8  5  2  3  1 
## 2 least connected regions:
## 15034 15061 with 1 link
## 1 most connected region:
## 15033 with 11 links
q1_existente <- nb2listw(queen.nb_existente)

class(queen.nb_existente) # Identificar el tipo de objeto
## [1] "nb"
class(q1_existente) # Identificar el tipo de objeto
## [1] "listw" "nb"
#k4.nb_existente= read.gwt2nb("./data/mx_city_metro_k4.gwt", region.id=muni$ID)
#summary(k4.nb_existente)

Crear pesos espaciales en R:

# queen weights:
queen.nb <- poly2nb(muni, row.names = muni$ID)
q1 <- nb2listw(queen.nb)
summary(queen.nb)
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 400 
## Percentage nonzero weights: 6.925208 
## Average number of links: 5.263158 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  2  3 11 13 16 12  8  5  2  3  1 
## 2 least connected regions:
## 36 59 with 1 link
## 1 most connected region:
## 35 with 11 links
summary(q1)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 400 
## Percentage nonzero weights: 6.925208 
## Average number of links: 5.263158 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  2  3 11 13 16 12  8  5  2  3  1 
## 2 least connected regions:
## 36 59 with 1 link
## 1 most connected region:
## 35 with 11 links
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 76 5776 76 32.83595 317.4801

nn = n*n S0 = suma global de los pesos

# rook weights:
rook.nb <- poly2nb(muni, row.names = muni$ID, queen=FALSE)
summary(rook.nb)
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 388 
## Percentage nonzero weights: 6.717452 
## Average number of links: 5.105263 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8 10 
##  2  4 10 17 15 10  7  7  4 
## 2 least connected regions:
## 36 59 with 1 link
## 4 most connected regions:
## 35 65 67 74 with 10 links

Mapear los pesos

plot(muni$geometry)
plot(queen.nb, st_geometry(st_centroid(muni)), cex = 0.6, pch = 20, col = 'red', add = TRUE)
plot(rook.nb, st_geometry(st_centroid(muni)), cex = 0.6, pch = 20, col = 'blue', add = TRUE)

nc.5nn<-knearneigh(st_geometry(st_centroid(muni)), k=5, longlat=TRUE)
nc.5nn.nb<-knn2nb(nc.5nn)
summary(nc.5nn.nb)
## Neighbour list object:
## Number of regions: 76 
## Number of nonzero links: 380 
## Percentage nonzero weights: 6.578947 
## Average number of links: 5 
## Non-symmetric neighbours list
## Link number distribution:
## 
##  5 
## 76 
## 76 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 with 5 links
## 76 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 with 5 links
knn5 <- nb2listw(nc.5nn.nb)
plot(muni$geometry)
xy <- st_coordinates(st_centroid(muni))[,1:2]
plot(nc.5nn.nb, xy, pch=20, col="blue", add=T)

Autocorrelación espacial en R

La dependencia espacial puede expresarse mediante autocorrelación espacial. La autocorrelación espacial revela valores similares en ubicaciones cercanas; se puede estudiar global o localmente. Concretamente, la autocorrelación espacial revelará si existe aleatoriedad espacial en la distribución de la inseguridad alimentaria entre municipios.

Para estudiar distintos patrones espaciales, la definición de pesos espaciales es crucial. La estructura de pesos espaciales determina la conectividad entre ubicaciones vecinas.

Estadísticas de autocorrelación espacial global

Ejemplos: Indice de Moran, Geary C, Getis-Ord, entre otras (mezclar similitud de atributos con pesos espaciales). \[ \sum_{ij}= w_{ij}f(x_ix_j) \]

Variable espacialmente rezagada Vale la pena examinar si hay autocorrelación espacial con variable de interés.

#Examinar si hay autocorrelación espacial con variable de interés
muni$Vx <- lag.listw(q1, muni$pqdwel_20)
plot(muni$pqdwel_20, muni$Vx)

Moran’s I (la estadística más usada): Los mapas muestran patrones de distribución espacial, dando indicios de presencia de autocorrelación espacial. Para formalmente evaluar la presencia de autocorrelación espacial global, utilizamos el índice univariado global de Moran.

\[ I= [\sum_i\sum_j w_{ij} z_iz_j/S_0]/[\sum_iz_i^2/N] \] donde \(S_0=\sum_i \sum_j w_{ij}\) y \(z_i = y_i-m_x\) es la desviación del promedio.

Es una estadístics de producto cruzado \((z_i z_j)\) similar a un coeficiente de correlación. En este caso, los valores dependen de una estructura de pesos \((w_{ij})\)

Hipótesis nula: La distribución de los datos es aleatoria en el espacio.

Consecuentemente, cuando la hipótesis nula es rechazada, podemos decir que hay diferencias espaciales en nuestros datos.

La autocorrelación positiva indica clustering, la autocorrelación negativa indica dispersión espaial.

Hacer estudios con distintos tipos de pesos para evaluar qué tanto se sostienen nuestras pruebas.

Estadística que puede rechazar la hipótesis nula por encontrar autocorrelación, no-normalidad. No sabemos bien qué es lo que está mal con nuestros datos.

Estadistica global de Moran I:

moran.test(muni$pqdwel_20, q1, alternative="two.sided")
## 
##  Moran I test under randomisation
## 
## data:  muni$pqdwel_20  
## weights: q1    
## 
## Moran I statistic standard deviate = 9.0233, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.644472659      -0.013333333       0.005314509

Prueba de Moran I:

moranTest <- moran.test(muni$pqdwel_20, q1)
moranTest
## 
##  Moran I test under randomisation
## 
## data:  muni$pqdwel_20  
## weights: q1    
## 
## Moran I statistic standard deviate = 9.0233, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.644472659      -0.013333333       0.005314509
moranTest$statistic
## Moran I statistic standard deviate 
##                           9.023318
moranTest$p.value
## [1] 9.123843e-20

Estadísticas de autocorrelación espacial local

Los clústeres se pueden identificar a través de indicadores locales de asociación espacial (LISA, Anselin 1995). Estos indicadores muestran las desviaciones de los patrones globales generales al mostrar la contribución de cada observación a la autocorrelación global general, lo que permite una visualización de los puntos calientes/fríos de la variable de interés Como se mencionó anteriormente, diferentes ponderaciones estudian la presencia de patrones mediante el uso de distintas técnicas de conglomerados, como el I de Moran local univariado.

En esta sección vamos a identificar valores atípicos (outliers) mediante las estadísticas de autocorrelación espacial local, identificando hot spots y cold spots.

Los ejemplos anteriores de estadísticas de autocorrelación espacial global están diseñados para encontrar si los datos (en su totalidad) tienen patrones de aleatoriedad espacial. Sin embargo, no muestran la ubicación de los clústers.

\[ \sum_{j}=w_{ij} f(x_ix_j)\]

Mapear la variable rezagada contra la variable

#Explicar cuadrantes
moran.plot(muni$pqdwel_20, q1)

Estimar autocorrelación local mediante la I de Moran univariada:

local <- localmoran(muni$pqdwel_20, listw = q1)

Mapear la estadistica local de Moran (Ii) - una LISA (Local Indicator of Spatial Autocorrelation):

moran.mapa <- cbind(muni, local)
tmap_mode("view")
tm_shape(moran.mapa) +
  tm_borders(alpha=0.7)+
  tm_fill(col = "Ii",
          alpha = 0.5,
          style = "quantile",
          title = "Estadistica Local de Moran I")

No olvidar que los aglomerados (clusters) que observamos incluyen a los vecinos.

Otra forma de estimar los cuadrantes de LISA para mapearlos y explorar nuestros datos:

#Examinemos aglomeraciones
quadrant <- vector(mode="numeric",length=nrow(local))

# Centra la variable de interés alrededor de promedio
m.pqdwel_20 <- muni$pqdwel_20 - mean(muni$pqdwel_20)     

# Centra estadística local de Moran's alrededor de promedio
m.local <- local[,1] - mean(local[,1])    

# Umbral de significancia
signif_1 <- 0.1 
signif_05 <- 0.05
signif_001 <- 0.001

# Construir cuadrantes para 0.1
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_1] <- 0   

# Armar mapa para 0.1
brks <- c(0,1,2,3,4)
colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")
plot(muni$geometry,border="lightgray",col=colors[findInterval(quadrant,brks,all.inside=FALSE)])
box() #Agregar caja
legend("bottomleft", legend = c("No significativo","Bajo-bajo","Bajo-alto","Alto-bajo","Alto-alto"),
       title="p-value = 0.1", fill=colors,bty="n")

# Forma alternativa
muni$LISA=quadrant
colores <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusteres <- c("No significativo","Bajo-bajo","Bajo-alto","Alto-bajo","Alto-alto")
tmap_mode("plot")
p1 = tm_shape(muni) +
  tm_fill(col = "LISA", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.1",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)
  

# Construir cuadrantes para 0.05
quadrant <- vector(mode="numeric",length=nrow(local))
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_05] <- 0  

muni$LISA=quadrant
p2 = tm_shape(muni) +
  tm_fill(col = "LISA", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.05",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

# Construir cuadrantes para 0.001
quadrant <- vector(mode="numeric",length=nrow(local))
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_001] <- 0  

muni$LISA=quadrant
p3 = tm_shape(muni) +
  tm_fill(col = "LISA", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.001",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

tmap_arrange(p1, p2, p3)

Ahora es importante estudiar los datos con otra definición de vecindad para comprobar la sensibilidad de nuestros análisis.

Con KNN = 5

Al considerar a los vecinos como todos los municipios abarcados dentro de una distancia determinada, los municipios más grandes tendrán menor conectividad y menos vecinos que los municipios más pequeños.

# con KNN = 5
local <- localmoran(muni$pqdwel_20, listw = knn5)

#Examinemos aglomeraciones
quadrant <- vector(mode="numeric",length=nrow(local))

# Centra la variable de interés alrededor de promedio
m.pqdwel_20 <- muni$pqdwel_20 - mean(muni$pqdwel_20)     

# Centra estadística local de Moran's alrededor de promedio
m.local <- local[,1] - mean(local[,1])    

# Umbral de significancia
signif_1 <- 0.1 
signif_05 <- 0.05
signif_001 <- 0.001

# Construir cuadrantes para 0.1
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_1] <- 0 

muni$LISA_kn=quadrant
colores <- c("#ffffff", "#2c7bb6", "#abd9e9", "#fdae61", "#d7191c")
clusteres <- c("No significativo","Bajo-bajo","Bajo-alto","Alto-bajo","Alto-alto")
tmap_mode("plot")
p1 = tm_shape(muni) +
  tm_fill(col = "LISA_kn", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.1",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)
  

# Construir cuadrantes para 0.05
quadrant <- vector(mode="numeric",length=nrow(local))
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_05] <- 0  

muni$LISA_kn=quadrant
p2 = tm_shape(muni) +
  tm_fill(col = "LISA_kn", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.05",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

# Construir cuadrantes para 0.001
quadrant <- vector(mode="numeric",length=nrow(local))
quadrant[m.pqdwel_20 >0 & m.local>0] <- 4  
quadrant[m.pqdwel_20 <0 & m.local<0] <- 1      
quadrant[m.pqdwel_20 <0 & m.local>0] <- 2
quadrant[m.pqdwel_20 >0 & m.local<0] <- 3
quadrant[local[,5]>signif_001] <- 0  

muni$LISA_kn=quadrant
p3 = tm_shape(muni) +
  tm_fill(col = "LISA_kn", style = "cat", palette = colores[c(sort(unique(quadrant)))+1],
          labels = clusteres[c(sort(unique(quadrant)))+1], popup.vars = c("ID")) +
  tm_borders(alpha=0.5)+
  tm_layout("LISA pvalue = 0.001",
            title.size = 1,
            legend.title.size = 1,
            legend.text.size = 0.6,
            legend.position = c("left","bottom"),
            legend.bg.color = "white",
            legend.bg.alpha = 1)

tmap_arrange(p1, p2, p3)

Heterogeneidad espacial en R

Se espera heterogeneidad espacial en México, donde algunas regiones son más heterogéneas que otras debido a las diferencias geográficas y económicas.

Imponer estructura de forma discreta (regímenes) o continua (smoothing, GWR - regresión geográficamente ponderada).

Creación de la variable dicotómica:

#Crear una variable dicotómica que sea 1 si el municipio es parte de la Ciudad de México y 0 si no.
muni$cdmx= ifelse(muni$state=="09", 1, 0)
municipios_4326$cdmx= ifelse(municipios_4326$state=="09", 1, 0)

# Ejemplo de manipulación de datos
#municipios_4326%>%group_by(cdmx)%>%summarise(n=n())

Mapear el régimen espacial creado:

ggplot(municipios_4326, aes(fill = cdmx)) + 
  geom_sf() +  # tells ggplot that geographic data are being plotted
  theme_minimal() + 
  labs(title = "Régimen espacial")

No sabemos por qué hay diferencias, pero resolvemos el problema de heterogeneidad.

Creación de variable dicotómica de densidad poblacional:

#Crear una variable dicotómica que sea 1 si el municipio es parte de la Ciudad de México y 0 si no.
municipios_4326$dens=as.numeric(municipios_4326$pop_20/st_area(municipios_4326)*1000)

mean(municipios_4326$dens)
## [1] 4.209716
municipios_4326$pop_dens= ifelse(municipios_4326$dens>4.222332, 1, 0)

#municipios_4326%>%group_by(pop)%>%summarise(n=n())
ggplot(municipios_4326, aes(fill = pop_dens)) + 
  geom_sf() +  # tells ggplot that geographic data are being plotted
  theme_minimal() + 
  labs(title = "Régimen espacial - densidad poblacional")