Directorio y librería

# Configuración
rm(list = ls())
# Paqueterías
pacman::p_load(tidyverse,
               ggplot2,
               viridis,
               sf,
               tmap,
               leaflet,
               #foreign,
               #expss,
               #fishualize,
               #raster,
               #cowplot
               #ggspatial,
               #colorspace
            )
## 
## Your package installed
## Warning in pacman::p_load(tidyverse, ggplot2, viridis, sf, tmap, leaflet, : Failed to install/load:

Shapes

Cargamos las Bases de datos geoespaciales (del INEGI) en R. Para más detalles de cómo funcionan los datos geoespaciales reomiendo el tutorial de Emmanuel Maruri

# Mapas a nivel Estado
mapa_mex_edos <- st_read("mg_sep2019_integrado/conjunto_de_datos/00ent.shp", 
                         options = "ENCODING=WINDOWS-1252")
## options:        ENCODING=WINDOWS-1252 
## Reading layer `00ent' from data source `/Users/Soporte/Desktop/Code/R Pro/Maps/mg_sep2019_integrado/conjunto_de_datos/00ent.shp' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 911292 ymin: 319149.1 xmax: 4082997 ymax: 2349615
## projected CRS:  MEXICO_ITRF_2008_LCC
glimpse(mapa_mex_edos)
## Rows: 32
## Columns: 4
## $ CVEGEO   <chr> "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", …
## $ CVE_ENT  <chr> "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", …
## $ NOMGEO   <chr> "Aguascalientes", "Baja California", "Baja California Sur", …
## $ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((2470518 115..., MULTIPOLYGON ((…
# Mapas a nivel Municipal
mapa_mex_mun <- st_read("mg_sep2019_integrado/conjunto_de_datos/00mun.shp", 
                        options = "ENCODING=WINDOWS-1252")
## options:        ENCODING=WINDOWS-1252 
## Reading layer `00mun' from data source `/Users/Soporte/Desktop/Code/R Pro/Maps/mg_sep2019_integrado/conjunto_de_datos/00mun.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2465 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 911292 ymin: 319149.1 xmax: 4082997 ymax: 2349615
## projected CRS:  MEXICO_ITRF_2008_LCC
glimpse(mapa_mex_mun)
## Rows: 2,465
## Columns: 5
## $ CVEGEO   <chr> "01001", "01002", "01003", "01004", "01005", "01006", "01007…
## $ CVE_ENT  <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", "01", …
## $ CVE_MUN  <chr> "001", "002", "003", "004", "005", "006", "007", "008", "009…
## $ NOMGEO   <chr> "Aguascalientes", "Asientos", "Calvillo", "Cosío", "Jesús Ma…
## $ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((2489073 111..., MULTIPOLYGON ((…

ENSANUT 2018

La Secretaría de Salud, el Instituto Nacional de Salud Pública (INSP) y el Instituto Nacional de Estadística y Geografía (INEGI), llevan a cabo el levantamiento de la Encuesta Nacional de Salud y Nutrición (ENSANUT) 2018, con el objeto de conocer el estado de salud y las condiciones nutricionales de la población en México.

A continuación, utilizamos los microdatos disponibles aquí para mapear la prevalencia de las principales enfermedades (diabetes, hipertensión y obesidad) en el país a nivel municipal.

En aspectos metodológicos la ENSANUT sigue la metodología de Unidad de Observación Vivienda, Hogares y el diseño estadístico de INEGI, por lo que podemos asociar los resultados presentados con los datos geoespaciales a nivel municipal. A continuación se muestra el proceso.

Datos de ensanut a nivel municipal

Importamos los microdatos del Componente de SALUD \(\Rightarrow\) Información sobre la vivienda a R desde el repositorio de Github donde están organizados.

ensanut_ap <- read.csv(
  file ="https://raw.githubusercontent.com/diego-eco/shiny_app_ensut/master/data/ensanut_areas_peq.csv", 
                       sep=",", 
                       colClasses=c(rep('character', 6), 'numeric','numeric','numeric')
                        )
names(ensanut_ap) <- c("CVEGEO", "ent", "ent_nom", "mun_clave", "mun_nom", 
                       "estimador", "Obesidad", "Hipertension", "Diabetes")
# Un vistazo rápido de la base
glimpse(ensanut_ap)
## Rows: 2,457
## Columns: 9
## $ CVEGEO       <chr> "01001", "01002", "01003", "01004", "01005", "01006", "0…
## $ ent          <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "…
## $ ent_nom      <chr> "Aguascalientes", "Aguascalientes", "Aguascalientes", "A…
## $ mun_clave    <chr> "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",…
## $ mun_nom      <chr> "Aguascalientes", "Asientos", "Calvillo", "Cosío", "Jesú…
## $ estimador    <chr> "Valor", "Valor", "Valor", "Valor", "Valor", "Valor", "V…
## $ Obesidad     <dbl> 31.5, 32.3, 40.0, 32.6, 34.7, 24.1, 31.7, 31.3, 32.0, 35…
## $ Hipertension <dbl> 14.9, 15.3, 13.8, 16.4, 12.4, 12.8, 14.3, 14.9, 14.7, 16…
## $ Diabetes     <dbl> 7.5, 8.0, 9.2, 7.4, 6.7, 7.4, 6.9, 6.2, 7.2, 7.6, 7.2, 9…

Notas que es importante que la columna de “mun” (municipios) del componente de salud de la ENSANUT la renombremos a “CVEGEO” para que coincida con la columna de los datos geoespaciales a nivel municipal.

Teniendo lo anterior listo, podemos unir las bases de datos para “añadir” los datos de ENSANUT a nuestra base de datos geoespaciales.

#Unir ambas bases
datos <- dplyr::left_join(mapa_mex_mun, ensanut_ap, by = c("CVEGEO"))
# Un vistazo rápido de la base
glimpse(datos)
## Rows: 2,465
## Columns: 13
## $ CVEGEO       <chr> "01001", "01002", "01003", "01004", "01005", "01006", "0…
## $ CVE_ENT      <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", "0…
## $ CVE_MUN      <chr> "001", "002", "003", "004", "005", "006", "007", "008", …
## $ NOMGEO       <chr> "Aguascalientes", "Asientos", "Calvillo", "Cosío", "Jesú…
## $ ent          <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "…
## $ ent_nom      <chr> "Aguascalientes", "Aguascalientes", "Aguascalientes", "A…
## $ mun_clave    <chr> "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",…
## $ mun_nom      <chr> "Aguascalientes", "Asientos", "Calvillo", "Cosío", "Jesú…
## $ estimador    <chr> "Valor", "Valor", "Valor", "Valor", "Valor", "Valor", "V…
## $ Obesidad     <dbl> 31.5, 32.3, 40.0, 32.6, 34.7, 24.1, 31.7, 31.3, 32.0, 35…
## $ Hipertension <dbl> 14.9, 15.3, 13.8, 16.4, 12.4, 12.8, 14.3, 14.9, 14.7, 16…
## $ Diabetes     <dbl> 7.5, 8.0, 9.2, 7.4, 6.7, 7.4, 6.9, 6.2, 7.2, 7.6, 7.2, 9…
## $ geometry     <MULTIPOLYGON [m]> MULTIPOLYGON (((2489073 111..., MULTIPOLYGO…

Podemos ver que en nuestra base “datos” tenemos listo los datos de diabetes, hipertensión y obesidad asociados a su variable “geometry” que nos permite mapear a nivel municipal.

Cuantiles

Antes de mapear, debemos tener en cuenta que tenemos el dato del porcentaje de población de 20 años y más con diagnóstico previp de la enfermedad para cada uno de los 2,465 municipios en México. No podemos mapear directamente estos datos porque necesitaríamos al menos 2,000 colores distintos en el mapa (asumiendo que algunos municipios tiene el mimso porcentaje que otros). Esto no es práctico en términos de cómputo ni en visualización de datos. Para más detalles sobre el propóstito de mapear ver ¿Por qué usar mapas?.

Por lo anterior, es importante antes de mapear dividir nuestras observaciones en cuantiles o rangos para poder asignar un color a cada uno. Esto se puede entender como “juntar” varias observaciones que compartan porcentajes cercanos, y asignarlas a un grupo/color.

Hay varias formas de agrupar los datos, a continuación mostramos un ejemplos con la variable de Diabetes:

Para extraer los cuantiles usamos las funciones: - pull() Parecido a select(), pero no igual - quantile()

Para crear las etiquetas utilizamos: - imap_chr() También vea aquí

# ¿Cuántas clases quiero?
no_classes <- 5
# Extraer cuantiles
cuantil_diab <- datos %>%
  pull(Diabetes) %>%
  quantile(probs = seq(0, 1, length.out = no_classes + 1),na.rm=T) %>%
  as.vector() # to remove names of quantiles, so idx below is numeric
# Así se crean las etiquetas
labels_diab <- imap_chr(cuantil_diab, function(., idx){
  return(paste0(round(cuantil_diab[idx] , 0),
                "%",
                " – ",
                round(cuantil_diab[idx + 1] , 0),
                "%"))
})
# Podemos observar los cuantiles que genera
cuantil_diab
## [1]  2.1  8.5  9.5 10.7 12.0 31.7
# Se elimina la última etiqueta del vector labels
# En caso contrario tendríamos una categoría que iría hasta NA
labels_diab <- labels_diab[1:length(labels_diab) - 1]
labels_diab
## [1] "2% – 8%"   "8% – 10%"  "10% – 11%" "11% – 12%" "12% – 32%"

Creamos la variable para cuantiles de Diabetes en la base de datos con la función cut() ver tamibén

# Crear la variable en datos
datos <- datos %>%
  mutate(q_Diabetes = cut(Diabetes,
                         breaks = cuantil_diab,
                         labels = labels_diab,
                         include.lowest = T))
glimpse(datos)
## Rows: 2,465
## Columns: 14
## $ CVEGEO       <chr> "01001", "01002", "01003", "01004", "01005", "01006", "0…
## $ CVE_ENT      <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", "0…
## $ CVE_MUN      <chr> "001", "002", "003", "004", "005", "006", "007", "008", …
## $ NOMGEO       <chr> "Aguascalientes", "Asientos", "Calvillo", "Cosío", "Jesú…
## $ ent          <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "…
## $ ent_nom      <chr> "Aguascalientes", "Aguascalientes", "Aguascalientes", "A…
## $ mun_clave    <chr> "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",…
## $ mun_nom      <chr> "Aguascalientes", "Asientos", "Calvillo", "Cosío", "Jesú…
## $ estimador    <chr> "Valor", "Valor", "Valor", "Valor", "Valor", "Valor", "V…
## $ Obesidad     <dbl> 31.5, 32.3, 40.0, 32.6, 34.7, 24.1, 31.7, 31.3, 32.0, 35…
## $ Hipertension <dbl> 14.9, 15.3, 13.8, 16.4, 12.4, 12.8, 14.3, 14.9, 14.7, 16…
## $ Diabetes     <dbl> 7.5, 8.0, 9.2, 7.4, 6.7, 7.4, 6.9, 6.2, 7.2, 7.6, 7.2, 9…
## $ geometry     <MULTIPOLYGON [m]> MULTIPOLYGON (((2489073 111..., MULTIPOLYGO…
## $ q_Diabetes   <fct> 2% – 8%, 2% – 8%, 8% – 10%, 2% – 8%, 2% – 8%, 2% – 8%, 2…

Ahora tenemos una nueva columna en neustra base de datos llamada q_Diabetes con la categoría asociada al municipio. Y nuestros 2,465 municipios ahora pueden agruparse en 5 clases y mapear de mejor forma, como ejemplo:

# Mapa prevalencia de diabetes con 5 clases a nivel municipal
 tm_shape(datos) +  
  tm_fill(col = "q_Diabetes")

En la siguiente sección presentamos los mapas con mayor detalle y estilo. Pero antes vamos a crear los cuantiles para nuestras otras dos variables: Hipertensión y Obesidad siguiendo el mismo procedimiento.

# ¿Cuántas clases quiero?
no_classes <- 5
# Extraer cuantiles
cuantil_hiper <- datos %>%
  pull(Hipertension) %>%
  quantile(probs = seq(0, 1, length.out = no_classes + 1),na.rm=T) %>%
  as.vector() # to remove names of quantiles, so idx below is numeric
# Así se crean las etiquetas
labels_hiper <- imap_chr(cuantil_hiper, function(., idx){
  return(paste0(round(cuantil_hiper[idx] , 0),
                "%",
                " – ",
                round(cuantil_hiper[idx + 1] , 0),
                "%"))
})

# Se elimina la última etiqueta
# En caso contrario sería hasta NA
labels_hiper <- labels_hiper[1:length(labels_hiper) - 1]

# Crear la variable en la base de datos
datos <- datos %>%
  mutate(q_Hipertension = cut(Hipertension,
                          breaks = cuantil_hiper,
                          labels = labels_hiper,
                          include.lowest = T))
# ¿Cuántas clases quiero?
no_classes <- 5
# Extraer cuantiles
cuantil_obes <- datos %>%
  pull(Obesidad) %>%
  quantile(probs = seq(0, 1, length.out = no_classes + 1),na.rm=T) %>%
  as.vector() # to remove names of quantiles, so idx below is numeric
# Así se crean las etiquetas
labels_obes <- imap_chr(cuantil_obes, function(., idx){
  return(paste0(round(cuantil_obes[idx] , 0),
                "%",
                " – ",
                round(cuantil_obes[idx + 1] , 0),
                "%"))
})

# Se elimina la última etiqueta
# En caso contrario sería hasta NA
labels_obes <- labels_obes[1:length(labels_obes) - 1]

# Crear la variable
datos <- datos %>%
  mutate(q_Obesidad = cut(Obesidad,
                              breaks = cuantil_obes,
                              labels = labels_obes,
                              include.lowest = T))
glimpse(datos)
## Rows: 2,465
## Columns: 16
## $ CVEGEO         <chr> "01001", "01002", "01003", "01004", "01005", "01006", …
## $ CVE_ENT        <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", …
## $ CVE_MUN        <chr> "001", "002", "003", "004", "005", "006", "007", "008"…
## $ NOMGEO         <chr> "Aguascalientes", "Asientos", "Calvillo", "Cosío", "Je…
## $ ent            <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1",…
## $ ent_nom        <chr> "Aguascalientes", "Aguascalientes", "Aguascalientes", …
## $ mun_clave      <chr> "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11…
## $ mun_nom        <chr> "Aguascalientes", "Asientos", "Calvillo", "Cosío", "Je…
## $ estimador      <chr> "Valor", "Valor", "Valor", "Valor", "Valor", "Valor", …
## $ Obesidad       <dbl> 31.5, 32.3, 40.0, 32.6, 34.7, 24.1, 31.7, 31.3, 32.0, …
## $ Hipertension   <dbl> 14.9, 15.3, 13.8, 16.4, 12.4, 12.8, 14.3, 14.9, 14.7, …
## $ Diabetes       <dbl> 7.5, 8.0, 9.2, 7.4, 6.7, 7.4, 6.9, 6.2, 7.2, 7.6, 7.2,…
## $ geometry       <MULTIPOLYGON [m]> MULTIPOLYGON (((2489073 111..., MULTIPOLY…
## $ q_Diabetes     <fct> 2% – 8%, 2% – 8%, 8% – 10%, 2% – 8%, 2% – 8%, 2% – 8%,…
## $ q_Hipertension <fct> 3% – 15%, 3% – 15%, 3% – 15%, 15% – 18%, 3% – 15%, 3% …
## $ q_Obesidad     <fct> 30% – 34%, 30% – 34%, 40% – 80%, 30% – 34%, 34% – 40%,…

Ahora tenemos listas nuestras 3 columnas en la base de datos llamadas q_Diabetes, q_Hipertension y q_Obesidad con la categoría asociada al municipio.

Notar que para las 3 enfermedades dividimos en 5 cateogrías al inicio de cada chunk, podemos repetir el proceso agreganod y quitando categorías o hacerlo interactivo en una Shiny app.

Mapas de prevalencia

Exiten muchas librerías para mapear en R, las principales son tmap y ggplot. Para entender tmap a detalle ver La función tmap en profundidad. A continuación mostramos los mapas de las 3 enfermedades con tmap.

Mapa de diabetes

# Mapa prevalencia de diabetes con 5 clases a nivel municipal
tm_shape(datos) +  
  tm_fill(col = "q_Diabetes",
          title = "% población > 20 años. Diagnóstico previo",
          palette = "Greens", 
          labels = labels_diab) + 
  tm_layout(main.title = "Prevalencia de diabetes en México", 
            main.title.position = "center",
            frame = F) +
  tm_credits("Fuente: Elaborado por Diego López con datos de la ENSANUT, INEGI.", 
             position = c("left", "bottom"))

Mapas de hipertensión

# Mapa prevalencia de hipertensión con 5 clases a nivel municipal
tm_shape(datos) +  
  tm_fill(col = "q_Hipertension",
          title = "% población > 20 años. Diagnóstico previo",
          palette = "Blues", 
          labels = labels_hiper) + 
  tm_layout(main.title = "Prevalencia de hipertensión en México", 
            main.title.position = "center",
            frame = F) +
  tm_credits("Fuente: Elaborado por Diego López con datos de la ENSANUT, INEGI.", 
             position = c("left", "bottom"))

Mapas con Obesidad

# Mapa prevalencia de Obesidad con 5 clases a nivel municipal
tm_shape(datos) +  
  tm_fill(col = "q_Obesidad",
          title = "% población > 20 años. Diagnóstico previo",
          palette = "Oranges", 
          labels = labels_obes) + 
  tm_layout(main.title = "Prevalencia de obesidad en México", 
            main.title.position = "center",
            frame = F) +
  tm_credits("Fuente: Elaborado por Diego López con datos de la ENSANUT, INEGI.", 
             position = c("left", "bottom"))

Color contínuo

A partir de aquí tomaremos la diabetes como nuestra enfermedad para generar los mapas.

Otro enfoque para mapear la prevalencia de enfermedades es utilizar un gradiente de color tomando las variable de prevalencia como un contínuo de porcentajes (entre 0-100).

Notar que para mapeo contínuo utilizamos directamente la variable diabetes y no q_Diabetes, porque ya no necesitamos los cuantiles generados. Podemos usar “order” o “cont” dependiendo de lo sesgado de la distribución con que estemos trabajando.

# Mapa prevalencia de diabetes con 5 clases a nivel municipal
tm_shape(datos) +  
  tm_fill(col = "Diabetes",
          style = "order",
          title = "% población > 20 años. Diagnóstico previo",
          palette = "-viridis") + 
  tm_layout(main.title = "Prevalencia de diabetes en México", 
            main.title.position = "center",
            frame = F) +
  tm_credits("Fuente: Elaborado por Diego López con datos de la ENSANUT, INEGI.", 
             position = c("left", "bottom"))

Filtrando para localidades

Podemos crear una base específica para una localidad, por ejemplo la CDMX con el número de entidad 9 (09, dependiendo de las etiquetas de ENSANUT o INEGI) y mapear la prevalencia de diabetes en las alcaldías.

# Filtramos CDMX
mapa_cdmx <- datos %>% filter(ent == "9")

tm_shape(mapa_cdmx) +  
  tm_fill(col = "q_Diabetes", 
          title = "% población > 20 años. Diagnóstico previo", 
          palette = "-RdYlGn", 
          labels = labels_diab) +
  tm_borders() + 
  tm_text(text = "NOMGEO",
          col = "black", 
          scale = 0.3, 
          remove.overlap = F,
          just = "top") +
  tm_layout(main.title = "Prevalencia de diabetes en CDMX", 
            main.title.position = "center", 
            main.title.size = 1.2,
            legend.outside = T, 
            frame = F) + 
  tm_credits("Fuente: Elaboración propia con \n datos de la ENSANUT, INEGI.",
             position = c("LEFT", "BOTTOM"))

Mapas interactivos

mapview

Podemos utilizar librerías como leaflet, mapview y plotly para generar mapas interactivos en HTML, para conocer más sobre estas herramientas ver Mapas interactivos rápidos con mapview ().

A continuación hacemos un mapa interactivo para la prevalencia de diabetes en el Valle de México.

mapa_valle_mexico <- datos %>% filter(ent %in% c("9", "15"))

library(mapview)
mapa_valle_mexico %>% mapview(color = "gray", lwd = 1, layer.name = "Municipios") +
 mapview(mapa_valle_mexico, zcol = "q_Diabetes")

leaflet

bins <- mapa_valle_mexico$q_Diabetes
pal <- colorBin("YlOrRd", domain = mapa_valle_mexico$q_Diabetes, bins = bins)


# Usamos pipes %>% para ir agregando las capas:
# leaflet(data = mapa_valle_mexico) %>% 
#  addProviderTiles(providers$CartoDB.Positron) %>% # Agregamos una capa de fondo de un provedor externo.
#  addCircles(col = ~pal(q_Diabetes), opacity = 0.9) %>% # Círculos con paleta según nbikes
#  addLegend(pal = pal, values = ~q_Diabetes) %>%  # Leyenda superior derecha usando la paleta pal
#  setView(lng = -0.1, 51.5, zoom = 12) %>%  # Vista inicial (posición y zoom)
#  addMiniMap() # Agrega el minimapa en la esquina inferior derecha (posición predeterminada

ggplot y plotly

Podemos crear mapas también a partir de la libería ggplot, ver a detalle en mapas en ggplot. Un beneficio de los mapas basados en ggplot2 es que se les puede dar fácilmente un nivel de interactividad cuando se imprimen usando la función ggplotly () del paquete plotly

map_valle <- ggplot(data = mapa_valle_mexico) +
            # Shape principal
            geom_sf(
                mapping = aes(
                    fill = q_Diabetes
                ),
                color = "white",
                size = 0
            ) +
            scale_fill_viridis(
                option = "cividis",
                name = " ",
                alpha = 0.8, 
                begin = 0.3, 
                end = 1,
                discrete = T, # True para clases discretas
                direction = -1, 
                guide = guide_legend(
                    keyheight = unit(5, units = "mm"),
                    title.position = "top",
                    reverse = T # El valor más alto hasta arriba
                )) +

            # Agregar títulos
            labs(x = NULL,
                 y = NULL,
                 title = "Prevalencia de Diabetes en el Valle de México",
                 subtitle = "% de población de 20 años y más con diagnóstico previo",
                 caption = "Fuente: Elaborado por Diego López con datos de INEGI/ENSANUT a nivel muncipal.",
                 fill = "%") +
             theme_bw()
plotly::ggplotly(map_valle)

  1. El Colegio de México, ↩︎