Situación problema

En México, se dice que los casos y muertes por COVID-19 sufridas durante la pandemia seguían ciertos patrones de comportamiento según el sector poblacional y sus características socio-económicas y geográficas. En este documento se pretende encontrar dichas características y validarlas con fundamentos estadísticos, para después poder identificar su comportamiento espacial, así como si tienen autocorrelación espacial. La base de datos que se utilizará con este propósito es la de ‘covid19_confirmados.csv’, usando como complemento la de ‘denue_hospitales.xlsx’.

El ESDA o Exploratory Spatial Data Analysis a diferencia del análisis de datos tradicional, busca encontrar una relación entre las variables de interés con ubicaciones. Nos sirve para poder identificar Autocorrelación Espacial, que confirma la presencia de variaciones espaciales. La autocorrelación espacial puede ser negativa o positiva, donde siendo positiva representa que las áreas cerca unas de las otras tienen valores similares, mientras que la negativa representa que las áreas de vecinos son diferentes.

Un ejemplo de autocorrelación espacial positiva son las características geográficas, es decir, un desierto en un territorio cuenta con la variable arena en toda su extensión territorial, por lo que se parecen entre más cerca están. Mientras tanto, un ejemplo de autocorrelación espacial negativa es

Carga de las bases de datos

Visualización de las bases de datos

str(hospitales)
## 'data.frame':    234303 obs. of  42 variables:
##  $ id        : int  9331895 6697315 6751037 536 8625042 16651 8463387 8984759 48213 8500272 ...
##  $ clee      : chr  "01001621111005611000000000U0" "01001621511001384002003549S9" "01001621511001441000003549S7" "01003624191000011001000000S3" ...
##  $ nom_estab : chr  "(H) CONSULTORIO AGUASCALIENTES, AGS." "176 SUCURSAL  EMILIANO ZAPATA SIGLO XXI" "184  SUCURSAL CONVENCI\xd3N ORIENTE  GRUPO DIAGN\xd3STICO M\xc9DICO PROA" "AA GRUPO CHICAGO 85" ...
##  $ raz_social: chr  "(H) CONSULTORIO AGUASCALIENTES, AGS. " "GRUPO DIAGNOSTICO MEDICO PROA SA DE CV" "GRUPO DIAGNOSTICO MEDICO PROA SA DE CV" "" ...
##  $ codigo_act: int  621111 621511 621511 624191 624191 624191 624211 623311 621331 621211 ...
##  $ nombre_act: chr  "Consultorios de medicina general del sector privado" "Laboratorios m\xe9dicos y de diagn\xf3stico del sector privado" "Laboratorios m\xe9dicos y de diagn\xf3stico del sector privado" "Agrupaciones de autoayuda para alcoh\xf3licos y personas con otras adicciones" ...
##  $ per_ocu   : chr  "0 a 5 personas" "31 a 50 personas" "0 a 5 personas" "11 a 30 personas" ...
##  $ tipo_vial : chr  "" "CALLE" "AVENIDA" "CALLE" ...
##  $ nom_vial  : chr  "" "GENERAL EMILIANO ZAPATA" "CONVENCI\xd3N DE 1914 ORIENTE" "ABASOLO" ...
##  $ tipo_v_e_1: chr  "" "CALLE" "AVENIDA" "CALLE" ...
##  $ nom_v_e_1 : chr  "" "JOS\xc9 GUADALUPE POSADA" "HALEY" "INDEPENDENCIA" ...
##  $ tipo_v_e_2: chr  "" "CALLE" "AVENIDA" "CALLE" ...
##  $ nom_v_e_2 : chr  "" "AZTECA" "PASEO DE LA CRUZ" "5 DE MAYO" ...
##  $ tipo_v_e_3: chr  "" "CALLE" "CALLE" "CALLE" ...
##  $ nom_v_e_3 : chr  "" "TALAMANTES" "V\xcdA L\xc1CTEA" "FRANCISCO JAVIER MINA" ...
##  $ numero_ext: int  NA 532 702 NA 201 400 202 111 45 1001 ...
##  $ letra_ext : chr  "" "" "" "SN" ...
##  $ edificio  : chr  "" "" "" "" ...
##  $ edificio_e: chr  "" "" "" "" ...
##  $ numero_int: int  NA 0 0 NA NA NA NA 0 NA 15 ...
##  $ letra_int : chr  "" "" "" "" ...
##  $ tipo_asent: chr  "" "COLONIA" "COLONIA" "COLONIA" ...
##  $ nomb_asent: chr  "" "SAN MARCOS" "JESUS GOMEZ PORTUGAL" "INDEPENDENCIA" ...
##  $ tipoCenCom: chr  "" "" "OTRO CONGLOMERADO" "" ...
##  $ nom_CenCom: chr  "" "" "LA LOMA" "" ...
##  $ num_local : chr  "" "" "702-A" "" ...
##  $ cod_postal: int  NA 20070 20250 20800 20326 20030 20140 20130 20674 20127 ...
##  $ cve_ent   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ entidad   : chr  "Aguascalientes" "Aguascalientes" "Aguascalientes" "Aguascalientes" ...
##  $ cve_mun   : int  1 1 1 3 1 1 1 1 6 1 ...
##  $ municipio : chr  "Aguascalientes" "Aguascalientes" "Aguascalientes" "Calvillo" ...
##  $ cve_loc   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ localidad : chr  "Aguascalientes" "Aguascalientes" "Aguascalientes" "Calvillo" ...
##  $ ageb      : chr  "731" "619" "2189" "187" ...
##  $ manzana   : int  3 9 5 24 15 8 4 17 13 19 ...
##  $ telefono  : chr  "" "" "" "4951077786" ...
##  $ correoelec: chr  "" "ERICK_ROQUE@PROA.COM.MX" "ERICK_ROQUE@PROA.COM.MX" "" ...
##  $ www       : chr  "" "WWW.CHOPO.COM.MX" "" "" ...
##  $ tipoUniEco: chr  "Fijo" "Fijo" "Fijo" "Fijo" ...
##  $ latitud   : num  21.9 21.9 21.9 21.8 21.9 ...
##  $ longitud  : num  -102 -102 -102 -103 -102 ...
##  $ fecha_alta: chr  "2020-11" "2016-01" "2016-10" "2010-07" ...
str(covid)
## 'data.frame':    2457 obs. of  42 variables:
##  $ cve_ent                       : int  1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 ...
##  $ mpio                          : chr  "Aguascalientes" "Asientos" "Calvillo" "Cosio" ...
##  $ poblacion_2022                : int  961977 50864 60760 16918 130184 50032 57981 9661 22743 21947 ...
##  $ hogrem2015                    : chr  "4.819" "12.867" "25.85" "14.294" ...
##  $ hogremjefmuj2015              : chr  "27.827" "22.054" "25.596" "21.171" ...
##  $ popnoafmed2015                : chr  "14.337" "5.774" "9.974" "5.283" ...
##  $ gini2015                      : chr  "0.392" "0.37" "0.375" "0.37" ...
##  $ popden2020                    : chr  "811.991" "94.087" "62.549" "131.885" ...
##  $ crimen_2018                   : num  7.36 7.81 8.4 17.86 5.14 ...
##  $ crimen_2019                   : num  8.15 5.78 8.31 0 10.98 ...
##  $ inclusion_fin_2019            : num  1.78 0 1.16 0 0.36 1.2 1.03 0 0 0 ...
##  $ porcentaje_pob_pobreza        : chr  "23.68" "40.13" "45.76" "37.04" ...
##  $ porcentaje_pob_pobreza_ext    : chr  "1.97" "4.14" "4.5" "3.38" ...
##  $ porcentaje_pob_servicios_salud: chr  "20.02" "16.47" "21.03" "17.59" ...
##  $ porcentaje_pob_acceso_ss      : chr  "37.77" "62.09" "76.15" "47.52" ...
##  $ pob_6.14_no_edu               : num  4.6 6.29 6.95 5.7 5.88 4.44 4.38 5.74 5.37 4.4 ...
##  $ rezago_social                 : num  -1.32 -0.86 -0.92 -1 -1.17 -1.17 -1.07 -0.96 -0.9 -0.83 ...
##  $ grado_rs                      : chr  "Muy bajo" "Muy bajo" "Muy bajo" "Muy bajo" ...
##  $ feb_2020                      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ march_2020                    : int  45 1 0 0 1 5 1 0 1 0 ...
##  $ april_2020                    : int  213 1 9 0 12 6 10 4 1 1 ...
##  $ may_2020                      : int  520 7 2 17 21 26 40 1 24 3 ...
##  $ june_2020                     : int  1168 21 25 25 42 57 88 18 43 3 ...
##  $ july_2020                     : int  1271 42 89 25 65 69 153 26 48 7 ...
##  $ august_2020                   : int  1385 17 18 17 71 40 85 6 24 4 ...
##  $ sept_2020                     : int  1364 13 41 2 65 82 47 3 14 4 ...
##  $ oct_2020                      : int  2388 46 64 6 93 109 74 9 19 6 ...
##  $ nov_2020                      : int  3965 94 185 13 104 192 73 13 24 2 ...
##  $ dic_2020                      : int  2370 73 116 2 73 87 16 14 18 11 ...
##  $ jan_2021                      : int  3271 64 186 9 78 90 42 9 22 7 ...
##  $ feb_2021                      : int  1498 42 151 7 73 46 34 7 8 4 ...
##  $ mar_2021                      : int  1310 21 73 0 38 34 10 4 8 5 ...
##  $ mar_2021.1                    : int  1310 21 73 0 38 34 10 4 8 5 ...
##  $ april_2021                    : int  826 14 19 2 27 15 11 5 2 6 ...
##  $ may_2021                      : int  401 6 3 1 14 3 13 0 0 3 ...
##  $ june_2021                     : int  221 2 0 4 8 3 3 1 1 0 ...
##  $ july_2021                     : int  1140 6 37 1 35 10 11 2 1 2 ...
##  $ august_2021                   : int  2944 49 82 4 116 74 52 12 19 3 ...
##  $ sept_2021                     : int  2121 23 22 1 95 38 79 11 34 11 ...
##  $ oct_2021                      : int  1560 33 57 4 78 39 66 5 12 22 ...
##  $ nov_2021                      : int  1126 16 32 4 45 38 61 2 10 14 ...
##  $ dic_2021                      : int  2067 14 68 4 95 18 30 0 6 11 ...

Uniendo las bases de datos.

#Preparing the variable 'clee' to match 'cve_ent' to do the joint
hospitales$cve_ent <- substr(hospitales$clee, 1, 5)
covid$cve_ent <- as.character(covid$cve_ent)
covid$cve_ent <- ifelse(str_length(covid$cve_ent) == 4, paste("0",covid$cve_ent, sep = ""),covid$cve_ent)
#Deleting all the health centers that aren't considered as hospitals or related to the treatment of Covid or related.
hospitalesFiltered <- filter(hospitales, codigo_act != 624191 & codigo_act !=
                               623311 & codigo_act !=
                               621331 & codigo_act !=
                               624198 & codigo_act !=
                               621311 & codigo_act !=
                               623991 & codigo_act !=
                               623992 & codigo_act !=
                               624221 & codigo_act !=
                               624112 & codigo_act !=
                               624111 & codigo_act !=
                               624311 & codigo_act !=
                               622311 & codigo_act !=
                               624411 & codigo_act !=
                               624412 & codigo_act !=
                               624122 & codigo_act !=
                               624222 & codigo_act !=
                               624312 & codigo_act !=
                               621411 & codigo_act !=
                               623312 & codigo_act !=
                               624121 & codigo_act !=
                               621412 & codigo_act !=
                               621312)
#Creating a Dataset that counts the number of hospitals in México
hospitalesMexico <- hospitalesFiltered %>% count(cve_ent,entidad)
hospitalesMexico <- hospitalesMexico[-c(2),]
#Merging the data bases into one 
covid$infectados2020 <- rowSums(covid[,c("feb_2020", "march_2020","april_2020","may_2020","june_2020","july_2020","august_2020","sept_2020","oct_2020","nov_2020","dic_2020")])
covid$infectados2021 <- rowSums(covid[,c("jan_2021", "feb_2021", "mar_2021","april_2021","may_2021","june_2021","july_2021","august_2021","sept_2021","oct_2021","nov_2021","dic_2021")])
joint <- merge(hospitales,covid, by.x = "cve_ent", by.y = "cve_ent", all.x = TRUE, all.y = FALSE)

Análisis exploratorio de la base de datos COVID

Buscando NA´s

Iniciamos con un análisis de NA´s que podrían estar presentes en la base de datos.

library(naniar)
gg_miss_var(covid, show_pct=TRUE)

No hay NA´s en la base de datos.

Transformación básica de la base de datos

Ahora, transformamos las columnas al tipo de dato que debe tener, además de crear dos nuevas columnas que tengan el total de casos de Covid por año. Así mismo, las columnas que sean representadas con porcentajes son divididas entre 100 para un mejor uso.

covid1 <-covid
covid1$casos2021 <- rowSums(covid[, 30:42])
covid1$casos2020 <- rowSums(covid[, 19:29])
covid1 <- covid1 %>% mutate_at(c("hogrem2015", "hogremjefmuj2015", "popnoafmed2015", "gini2015", "porcentaje_pob_pobreza", "porcentaje_pob_pobreza_ext", "porcentaje_pob_servicios_salud", "porcentaje_pob_acceso_ss", "popden2020"), as.numeric)
covid1 <- na.omit(covid1)
covid1$porcentaje_pob_pobreza <- covid1$porcentaje_pob_pobreza/100
covid1$porcentaje_pob_acceso_ss <- covid1$porcentaje_pob_acceso_ss/100
covid1$porcentaje_pob_pobreza_ext <- covid1$porcentaje_pob_pobreza_ext/100
covid1$porcentaje_pob_servicios_salud <- covid1$porcentaje_pob_servicios_salud/100

Gráfica de correlación entre variables en la base de datos.

library(dlookr)
corr_covid <- covid1 %>% select(hogrem2015, hogremjefmuj2015, popnoafmed2015, gini2015, porcentaje_pob_pobreza, porcentaje_pob_pobreza_ext, porcentaje_pob_servicios_salud, porcentaje_pob_acceso_ss, popden2020, casos2020)
correlate(corr_covid, hogrem2015, hogremjefmuj2015, popnoafmed2015, gini2015, porcentaje_pob_pobreza, porcentaje_pob_pobreza_ext, porcentaje_pob_servicios_salud, porcentaje_pob_acceso_ss, popden2020, casos2020) %>%  plot()

Como se puede apreciar en la anterior gráfica de correlación, las columnas casos2020, porcentaje_pob_pobreza y porcentaje_pob_acceso_ss, son las columnas que más correlación presentan, además de la pobreza extrema y densidad de población, por lo que estás serán las analizadas con más detalle.

Variables seleccionadas

  • Municipio se refiere ala división territorial administrativa que está regida por un ayuntamiento, por estado.

  • Popden2020 se refiere a la concentración de población (densidad) en una locación geográfica; se utiliza para cuantificar información demográfica en relación a ecosistemas, salud, infraestructura y distintos patrones.

  • Porcentaje_pob_pobreza explica en porcentaje, la pobreza de la población en las áreas especificadas de México. Según el Consejo Nacional de Evaluación de la Política de Desarrollo Social (CONEVAL), la pobreza por ingresos se define por tres líneas. La línea de pobreza alimentaria habla de la capacidad para obtener una canasta básica alimentaria, aquella de capacidades a la suficiencia del ingreso disponible para adquirir el valor de la canasta alimentaria y efectuar los gastos necesarios en salud y educación y finalmente la de patrimonio que habla de la insuficiencia de éstos dos así como de vivienda, transporte y bienes y servicios.

  • Porcentaje_pob_pobreza_ext se refiere a la situación de pobreza extrema al tener tres o más carencias (de seis) del Índice de Privación Social y estar por debajo de la línea de bienestar (mínima).

  • Porcentaje_pob_servicios_salud se refiere al porcentaje de la población por área, que tiene acesso a servicios de salud pública y privada

  • Porcentaje_pob_acceso_ss se refiere al porcentaje de la población que tiene acceso al seguro social.

Así mismo, se utilizarán las variables HospitalesHabitantes y casos2020, las cuales serán calculadas más adelante.

  • HospitalesHabitantes: Indica el número de hospitales cada 10 mil habitantes.

  • Casos2020: Indica el número de casos de COVID 19 en el año 2020.

Análisis de las variables seleccionadas

Gráficas de Normalidad

plot_normality(covid1, casos2020, porcentaje_pob_pobreza, porcentaje_pob_acceso_ss)

Para poder identificar la distribución de los datos en las variables de interés se hacen estás 12 gráficas, dando como resultado de solo la variable de casos2020 seguiría una distribución normal solo si se le aplica una transformación logarítmica.

Análisis Bivariado entre Población con Pobreza y Casos2020

target_var_1<-target_by(covid1, porcentaje_pob_pobreza)
relationship_1<-relate(target_var_1,  casos2020)
plot(relationship_1)

Análisis Bivariado entre Accesso a SS y Casos2020

target_var_2<-target_by(covid1, porcentaje_pob_acceso_ss)
relationship_2<-relate(target_var_2, casos2020)
plot(relationship_2)

Gráfica: Relación entre pob_pobreza y accesso_ss en casos2020

ggplot(covid1, aes(x=porcentaje_pob_pobreza, y=porcentaje_pob_acceso_ss, size=casos2020)) +
  geom_point() +
  scale_x_log10()

Como se puede ver en esta última gráfica de exploración de los datos, entre menor es el acceso al Seguro Social y mayor la pobreza, incremente considerablemente la cantidad de casos de Covid, confirmando así lo que hemos visto en todas nuestras gráficas. Por lo tanto, las variables porcentaje_pob_pobreza y porcentaje_pob_pobreza son explicativas y nos pueden servir más adelante cuando veamos la distribución de hospitales en México.

Visualizando la cantidad de hospitales en las regiones de México.

#Creating subsets for the regions in México
hospitalesCentroNorte <- filter(hospitalesFiltered, entidad == "Aguascalientes"
                                          | entidad == "Baja California Sur"
                                          | entidad == "Colima"
                                          | entidad == "Durango"
                                          | entidad == "Jalisco"
                                          | entidad == "Michoac\xe1n de Ocampo"
                                          | entidad == "Nayarit"
                                          | entidad == "San Luis Potos\xed"
                                          | entidad == "Sinaloa"
                                          | entidad == "Zacatecas")

hospitalesCentro <- filter(hospitalesFiltered, entidad == "Ciudad de M\xe9xico"
                                    | entidad == "M\xe9xico"
                                    | entidad == "Guanajuato"
                                    | entidad == "Hidalgo"
                                    | entidad == "Morelos"
                                    | entidad == "Puebla"
                                    | entidad == "Quer\xe9taro"
                                    | entidad == "Tlaxcala")

hospitalesNorte <- filter(hospitalesFiltered, entidad == "Baja California" 
                          | entidad == "Chihuahua"
                          | entidad == "Coahuila de Zaragoza"
                          | entidad == "Nuevo Le\xf3n"
                          | entidad == "Sonora"
                          | entidad == "Tamaulipas")

hospitalesSur <- filter(hospitalesFiltered, entidad == "Campeche"
                                  | entidad == "Chiapas"
                                  | entidad == "Guerrero"
                                  | entidad == "Oaxaca"
                                  | entidad == "Quintana Roo"
                                  | entidad == "Tabasco"
                                  | entidad == "Veracruz de Ignacio de la Llave"
                                  | entidad == "Yucat\xe1n")
#Adding the amount of hospitals for region
jointCentroNorte <- merge(hospitalesCentroNorte,covid, by.x = "cve_ent", by.y = "cve_ent", all.x = TRUE, all.y = FALSE)
jointCentro <- merge(hospitalesCentro,covid, by.x = "cve_ent", by.y = "cve_ent", all.x = TRUE, all.y = FALSE)
jointNorte <- merge(hospitalesNorte,covid, by.x = "cve_ent", by.y = "cve_ent", all.x = TRUE, all.y = FALSE)
jointSur <- merge(hospitalesSur,covid, by.x = "cve_ent", by.y = "cve_ent", all.x = TRUE, all.y = FALSE)

Mapa: Cantidad de hospitales en México

qmplot(longitud, latitud, data = joint, colour = I('red'), size = I(0.8), darken = .3)

Mapa: Cantidad de hospitales en el Centro Norte de México

qmplot(longitud, latitud, data = jointCentroNorte, colour = I('red'), size = I(0.8), darken = .3)

Mapa: Cantidad de hospiales en el Centro de México

qmplot(longitud, latitud, data = jointCentro, colour = I('red'), size = I(0.8), darken = .3)

Mapa: Cantidad de hospitales en el Sur de México

qmplot(longitud, latitud, data = jointSur, colour = I('red'), size = I(0.8), darken = .3)

Mapa: Cantidad de hospitales en el Norte de México

qmplot(longitud, latitud, data = jointNorte, colour = I('red'), size = I(0.8), darken = .3)

densidad <- readxl::read_excel("C:\\Users\\art191127\\Downloads\\Población_07.xlsx")
superficie <- readxl::read_excel("C:\\Users\\art191127\\Downloads\\inegi_superficie_territorial.xlsx",sheet = "Distribución territorial")

Creando la variable Hospitales x Habitantes

superficie <- na.omit(superficie)
colnames(superficie) <- c("Entidad federativa", "Municipio", "Superficie (Km2)", "Porcentaje de la superficie estatal","Demsidad de población (hab./km2)","Total de municipios o demarcaciones territoriales / localidades")
superficie$`Entidad federativa` <- substr(superficie$`Entidad federativa`,1,2)
superficie$Municipio <- substr(superficie$Municipio,1,3)

superficie$cve_ent <- paste(superficie$`Entidad federativa`,superficie$Municipio,sep = "")

#Adding population density into hospitalesMexico dataset

hospitalesMexico <- merge(hospitalesMexico,superficie,by.x = "cve_ent", by.y = "cve_ent")
hospitalesMexico <- merge(hospitalesMexico, covid %>% select(cve_ent, poblacion_2022), by.x = "cve_ent", by.y = "cve_ent")

#Deleting extra municipalities 

hospitalesMexico <- hospitalesMexico[-c(2457, 2458), ]

#Creating the columns for KM2 extension and population density

hospitalesMexico$hospitalesKm2 <- hospitalesMexico$n/hospitalesMexico$`Superficie (Km2)`
hospitalesMexico$hospitalespHabitantes <- (hospitalesMexico$n/hospitalesMexico$poblacion_2022)*10000
map <- st_read("C:\\Users\\art191127\\Documents\\Tec de Monterrey\\Semestre 8\\Clase 2\\Saucedo\\shp_mx_mpios\\mx_mpios.shp")
## Reading layer `mx_mpios' from data source 
##   `C:\Users\art191127\Documents\Tec de Monterrey\Semestre 8\Clase 2\Saucedo\shp_mx_mpios\mx_mpios.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
map_deprctd <- readShapePoly("C:\\Users\\art191127\\Documents\\Tec de Monterrey\\Semestre 8\\Clase 2\\Saucedo\\shp_mx_mpios\\mx_mpios.shp", IDvar = "IDUNICO", proj4string = CRS("+proj=longlat"), sf_use_s2(FALSE))
## Shapefile type: Polygon, (5), # of Shapes: 2456
hospitalesMexico$cve_ent <- as.numeric(hospitalesMexico$cve_ent)
merged_data <- merge(map,hospitalesMexico, by.x = "IDUNICO", by.y = "cve_ent", all.y = TRUE)

Explorando las variables seleccionadas en México

Mapa: Cantidad de Hospitales cada 10 mil habitantes.

municipiosHH <- ggplot() + 
  geom_sf(data = merged_data, aes(fill = hospitalespHabitantes)) +
  scale_fill_gradient(low = "#fff7ec", high = "#7F0000", limits = c(0.000009, 28))

municipiosHH

Iniciando por la cantidad de hospitales cada 10 mil habitantes, podemos observar como en el norte predomina una mayor cantidad de hospitales para sus habitantes, mientras que en el sur, hay un concentración mucho menor.

# Transforming selected variables to plot them 
covid <- covid %>% mutate_at(c("hogrem2015", "hogremjefmuj2015", "popnoafmed2015", "gini2015", "porcentaje_pob_pobreza", "porcentaje_pob_pobreza_ext", "porcentaje_pob_servicios_salud", "porcentaje_pob_acceso_ss", "popden2020"), as.numeric)
#covid1 <- na.omit(covid1)
# Dividing the percentage variables between 100 for a better interpretation
covid$porcentaje_pob_pobreza <- covid$porcentaje_pob_pobreza/100
covid$porcentaje_pob_acceso_ss <- covid$porcentaje_pob_acceso_ss/100
covid$porcentaje_pob_pobreza_ext <- covid$porcentaje_pob_pobreza_ext/100
covid$porcentaje_pob_servicios_salud <- covid$porcentaje_pob_servicios_salud/100
covid$cve_ent <- as.character(covid$cve_ent)
colnames(map)[3] ="cve_ent"
map$cve_ent <- as.character(map$cve_ent)
# Deleting the extra row
map <- map[-2456,]
map$cve_ent <- as.character(map$cve_ent)
# Matching the key of the dataset 'map' with the key of 'covid'
map$cve_ent <- ifelse(str_length(map$cve_ent) == 4, paste("0",map$cve_ent, sep = ""),map$cve_ent)
connectivityM <- merge(map, covid, by = "cve_ent")
# Selecting only the interest variables 
vector <- as.vector(c("cve_ent","porcentaje_pob_pobreza_ext", "porcentaje_pob_pobreza", "porcentaje_pob_servicios_salud", "rezago_social", "popden2020","infectados2020", "crimen_2019", "geometry"))
#vector <- as.vector(c("porcentaje_pob_pobreza", "rezago_social", "feb_2020", "march_2020","april_2020","may_2020","june_2020","july_2020","august_2020","sept_2020","oct_2020","nov_2020","dic_2020", "geometry"))
connectivityM <- connectivityM[,vector]  
#mex_cov <- function(df, fill) {
  #ggplot() + 
  #geom_sf(data = df, aes(fill = fill)) +
  #scale_fill_gradient(low = "#fff7ec", high = "#7F0000", ggtitle("Porcentaje de Pobreza Ext en México"))
#}

#mex_cov(connectivityM, porcentaje_pob_pobreza_ext)

####
municipiosPX <- ggplot() + 
  geom_sf(data = connectivityM, aes(fill = porcentaje_pob_pobreza_ext)) +
  scale_fill_gradient(low = "#fff7ec", high = "#7F0000", ggtitle("Porcentaje de Pobreza Ext en México"))

municipiosPX

Mapa: Pobreza extrema en México

municipiosPX <- ggplot() + 
  geom_sf(data = connectivityM, aes(fill = porcentaje_pob_pobreza_ext)) +
  scale_fill_gradient(low = "#fff7ec", high = "#7F0000", ggtitle("Porcentaje de Pobreza Ext en México"))

municipiosPX

En este mapa se puede observar un tendencia contraria al anterior, donde la pobreza extrema está concentrada más en el centro y sur del país.

Mapa: Porcentaje de Pobreza en México

municipiosPX <- ggplot() + 
  geom_sf(data = connectivityM, aes(fill = porcentaje_pob_pobreza)) +
  scale_fill_gradient(low = "#fff7ec", high = "#7F0000", ggtitle("Porcentaje de Pobreza en México"))

municipiosPX

Al igual que el anterior, en este caso podemos observar que la pobreza es más común en el sur del país, aunque a diferencia de la pobreza extrema, en este caso podemos ver más casos en el norte.

Mapa: Porcentaje de Carencia a Servicios de Salud en México

municipiosPX <- ggplot() + 
  geom_sf(data = connectivityM, aes(fill = porcentaje_pob_servicios_salud)) +
  scale_fill_gradient(low = "#fff7ec", high = "#7F0000", ggtitle("Por de Carencia a SS en México"))

municipiosPX

Siguiendo con la tendencia, el sur nuevamente vuelve a ser el que más problemas presenta, ahora en términos de carencia en Servicios de Salud

Mapa: Rezago social en México

municipiosPX <- ggplot() + 
  geom_sf(data = connectivityM, aes(fill = rezago_social)) +
  scale_fill_gradient(low = "#fff7ec", high = "#7F0000", ggtitle("Rezago Social en México"))

municipiosPX

En el caso del rezago social, es más difícil identificar una tendencia, sin embargo, hay unos pocos municipios que se destacan claramente.

Mapa: Densidad de Población en México

Al hacer este mapa salió a la luz que la variable de densidad de población presentaba problemas, por lo que se hicieron 2 boxplots previas.

boxplot(connectivityM$popden2020)

En la primera podemos apreciar como esta variable cuenta con outliers muy altos, probablemente el más alto siendo el Estado de México, sin embargo, al no poder eliminar este municipio, se optó por hacer una transformación logarítmica.

boxplot(log(connectivityM$popden2020))

Como se puede ver después de la transformación, la variable cuenta con una distribución mucho más sana.

municipiosPX <- ggplot() + 
  geom_sf(data = connectivityM, aes(fill = log(popden2020))) +
  scale_fill_gradient(low = "#fff7ec", high = "#7F0000", ggtitle("Densidad de Población"))

municipiosPX

Aunque no se aprecia en su totalidad, el centro de México es el más poblado.

Mapa: Crimen en México

Esta variable también presento un problema al momento de ser graficada, por lo que se hizo un summary de la variable.

summary(connectivityM$crimen_2019)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00   11.49   20.29   26.92  551.82

Como se puede ver, la variable cuenta con el valor 0, el cual es común en varios municipios, por lo que es dificil saber si se trata de un error o es intencional. En cualquier caso se graficó con los ceros.

municipiosPX <- ggplot() + 
  geom_sf(data = connectivityM, aes(fill = crimen_2019)) +
  scale_fill_gradient(low = "#fff7ec", high = "#7F0000", limits = c(0, 27), ggtitle("Crimen en México"))

municipiosPX

Dado la naturaleza de la variable, el mapa presenta muchas zonas grises que representan los ceros en la variable, por lo que al haber tantos, la interpretación se vuelve confusa, por lo que está variable no seguirá siendo tomada en cuenta.

SPATIAL WEIGHT MATRIX MÉXICO

Para iniciar con nuestro análisis de Autocorrelación Espacial, decidimos hacer 6 gráficas de Spatial Weight Matrix, donde se busca identificar el comportamiento espacial de los datos, en este caso en tipo ‘rook’ lo que permite identificar vecinos incluyendo los vértices comunes. Por lo tanto observamos 6 variables graficadas, porcentaje_pob_pobreza_ext, porcentaje_pob_pobreza, porcentaje_pob_servicios_salud, rezago_social, popden2020 y crimen_2019.

hospitalesMexico$cve_ent <- as.character(hospitalesMexico$cve_ent)
hospitalesMexico$cve_ent <- ifelse(str_length(hospitalesMexico$cve_ent) == 4, paste("0",hospitalesMexico$cve_ent, sep = ""),hospitalesMexico$cve_ent)
connectivityM <- merge(connectivityM, hospitalesMexico, by = "cve_ent")
vector <- as.vector(c("porcentaje_pob_pobreza_ext", "porcentaje_pob_pobreza", "porcentaje_pob_servicios_salud", "rezago_social", "popden2020","hospitalespHabitantes", "infectados2020", "geometry"))
connectivityM <- connectivityM[,vector]  
#map_deprctd <- map_deprctd[-2457,]
connectivityM <- na.omit(connectivityM)

Mapas: SWM México

lmat_c<-coordinates(map_deprctd)
map.centroid_c<-coordinates(map_deprctd)  
connectivityM.link_a_rook<-poly2nb(connectivityM,queen=T)
connectivityM.linkW_a_rook<-nb2listw(connectivityM.link_a_rook, style="W")
plot(connectivityM,border="white",axes=TRUE,las=1)
plot(connectivityM,col="grey",border=grey(0.11),axes=T,add=T) 
plot(connectivityM.linkW_a_rook,coords=map.centroid_c,pch=19,cex=0.1,col="red",add=T)

Al igual que con los mapas de exploración, podemos observar como las variables de pobreza, acceso a servicios de salud y rezago social presentan los

hospitalesMexico_centroid<-coordinates(map_deprctd)
map_nb<-poly2nb(map_deprctd)
hospitalesMexico_map.linkW<-nb2listw(map_nb, style="W")   
plot(map_deprctd,border="blue",axes=FALSE,las=1, main="MX Spatial Connectivity Matrix")
plot(map_deprctd,col="grey",border=grey(0.9),axes=T,add=T) 
plot(hospitalesMexico_map.linkW,coords=hospitalesMexico_centroid,pch=19,cex=0.1,col="red",add=T)

Autocorrelación Global para las variables

Moran test: Porcentaje de Pobreza Extrema

library(spdep)

moran.test(connectivityM$porcentaje_pob_pobreza_ext, nb2listw(connectivityM.link_a_rook))
## 
##  Moran I test under randomisation
## 
## data:  connectivityM$porcentaje_pob_pobreza_ext  
## weights: nb2listw(connectivityM.link_a_rook)    
## 
## Moran I statistic standard deviate = 56.891, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.6927728299     -0.0004078303      0.0001484580
moran.plot(connectivityM$porcentaje_pob_pobreza_ext, nb2listw(connectivityM.link_a_rook))

Tanto la gráfica como el resultado del Moran Test indican autocorrelación espacial positiva, dado que el Moran statistic es de 0.69 y el p-value es significativo.

Además, dado que el Moran I statistic es mayor a 0.5, se puede inferir que existen clusters a nivel local que concentran mayores cantidades de Pobreza Extrema, a comparación de lo que se podría observar en otros lugares del país.

Moran test: Porcentaje de Pobreza

moran.test(connectivityM$porcentaje_pob_pobreza, nb2listw(connectivityM.link_a_rook))
## 
##  Moran I test under randomisation
## 
## data:  connectivityM$porcentaje_pob_pobreza  
## weights: nb2listw(connectivityM.link_a_rook)    
## 
## Moran I statistic standard deviate = 62.273, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.7587271433     -0.0004078303      0.0001486072
moran.plot(connectivityM$porcentaje_pob_pobreza, nb2listw(connectivityM.link_a_rook))

Tanto la gráfica como el resultado del Moran Test indican autocorrelación espacial positiva, dado que el Moran statistic es de 0.75 y el p-value es significativo.

De igual forma que el anterior el Moran I statistic es mayor a 0.5, por lo que es un valor alto y se puede inferir que existen clusters locales que tienen mayores niveles de pobreza a comparación de otros lugares en el país.

Moran test: Porcentaje de Acceso a Servicios de Salud

moran.test(connectivityM$porcentaje_pob_servicios_salud, nb2listw(connectivityM.link_a_rook))
## 
##  Moran I test under randomisation
## 
## data:  connectivityM$porcentaje_pob_servicios_salud  
## weights: nb2listw(connectivityM.link_a_rook)    
## 
## Moran I statistic standard deviate = 34.389, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.4187152972     -0.0004078303      0.0001485366
moran.plot(connectivityM$porcentaje_pob_servicios_salud, nb2listw(connectivityM.link_a_rook))

Tanto la gráfica como el resultado del Moran Test indican autocorrelación espacial positiva, dado que el Moran statistic es de 0.42 y el p-value es significativo.

Moran I statistic menor a 0.5, no se infiere la existencia de clusters locales.

Moran test: Porcentaje de Rezago Social

moran.test(connectivityM$rezago_social, nb2listw(connectivityM.link_a_rook))
## 
##  Moran I test under randomisation
## 
## data:  connectivityM$rezago_social  
## weights: nb2listw(connectivityM.link_a_rook)    
## 
## Moran I statistic standard deviate = 56.589, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.6889664298     -0.0004078303      0.0001484024
moran.plot(connectivityM$rezago_social, nb2listw(connectivityM.link_a_rook))

Tanto la gráfica como el resultado del Moran Test indican autocorrelación espacial positiva, dado que el Moran statistic es de 0.68 y el p-value es significativo.

Moran I statistic mayor a 0.5, por lo que se infiere la existencia de clusters locales que concentran mayor cantidades de Rezago Social a comparación de otros lugares del país.

Moran test: Densidad de Población

moran.test(log(connectivityM$popden2020), nb2listw(connectivityM.link_a_rook))
## 
##  Moran I test under randomisation
## 
## data:  log(connectivityM$popden2020)  
## weights: nb2listw(connectivityM.link_a_rook)    
## 
## Moran I statistic standard deviate = 60.396, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.7356359156     -0.0004078303      0.0001485205
moran.plot(log(connectivityM$popden2020), nb2listw(connectivityM.link_a_rook))

Tanto la gráfica como el resultado del Moran Test indican autocorrelación espacial positiva, dado que el Moran statistic es de 0.73 y el p-value es significativo.

Moran I statistic mayor a 0.5, por lo que se infiere la existencia de clusters locales que concentran mayor cantidades de Densidad Poblacional a comparación de otros lugares en el país.

Moran test: Hospitales por Habitantes

moran.test(connectivityM$hospitalespHabitantes, nb2listw(connectivityM.link_a_rook))
## 
##  Moran I test under randomisation
## 
## data:  connectivityM$hospitalespHabitantes  
## weights: nb2listw(connectivityM.link_a_rook)    
## 
## Moran I statistic standard deviate = 23.349, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2822790706     -0.0004078303      0.0001465748
moran.plot(connectivityM$hospitalespHabitantes, nb2listw(connectivityM.link_a_rook))

Tanto la gráfica como el resultado del Moran Test indican autocorrelación espacial positiva, dado que el Moran statistic es de 0.28 y el p-value es significativo.

Moran I statistic menor a 0.5, por lo que se infiere que no existen clusters locales.

Local Cluster Analysis

Tomando en consideración los hallazgos previos, hacemos un análisis de clusters locales LISA para las variables seleccionadas.

Mapa: LISA Analysis para Pobreza Extrema

map <- map[-2455,]
map <- map[-2454,]
queen_w<-queen_weights(map)
connectivityM1 <- connectivityM

lisa_pppe<-local_moran(queen_w, connectivityM["porcentaje_pob_pobreza_ext"]) 

connectivityM1$porcentaje_pob_pobreza_ext<-as.factor(lisa_pppe$GetClusterIndicators())
levels(connectivityM1$porcentaje_pob_pobreza_ext)<-lisa_pppe$GetLabels() 

ggplot(data=connectivityM1) + 
  geom_sf(aes(fill=porcentaje_pob_pobreza_ext)) + 
  ggtitle(label="Porcentaje de Pobreza Extrema")

Como se puede ver en el mapa, no hay muchos clusters autocorrelación espacial positiva o de cualquier tipo en general.

Mapa: LISA Analysis para Pobreza

lisa_pppe<-local_moran(queen_w, connectivityM["porcentaje_pob_pobreza"]) 

connectivityM1$porcentaje_pob_pobreza<-as.factor(lisa_pppe$GetClusterIndicators())
levels(connectivityM1$porcentaje_pob_pobreza)<-lisa_pppe$GetLabels() 

ggplot(data=connectivityM1) + 
  geom_sf(aes(fill=porcentaje_pob_pobreza)) + 
  ggtitle(label="Porcentaje de Pobreza")

En este caso, en la zona norte del país se puede ver que en varios municipios, si la pobreza el baja, a su alrededor también tiende a serlo.

Mapa: LISA Analysis para Acceso a Servicios de Salud

lisa_pppe<-local_moran(queen_w, connectivityM["porcentaje_pob_servicios_salud"]) 

connectivityM1$porcentaje_pob_servicios_salud<-as.factor(lisa_pppe$GetClusterIndicators())
levels(connectivityM1$porcentaje_pob_servicios_salud)<-lisa_pppe$GetLabels() 

ggplot(data=connectivityM1) + 
  geom_sf(aes(fill=porcentaje_pob_servicios_salud)) + 
  ggtitle(label="Porcentaje Acceso a Servicios de Salud")

Sin muchos patrones, más que en el norte podemos ver pocos municipios con LL que significa que si hay baja cantidad de falta a acceso a servicios de salud, tiende a también haberla en la cercanía. Pero como se esperaba desde el Moran Test, realmente no existen clusters.

Mapa: LISA Analysis para Rezago Social

lisa_pppe<-local_moran(queen_w, connectivityM["rezago_social"]) 

connectivityM1$rezago_social<-as.factor(lisa_pppe$GetClusterIndicators())
levels(connectivityM1$rezago_social)<-lisa_pppe$GetLabels() 

ggplot(data=connectivityM1) + 
  geom_sf(aes(fill=rezago_social)) + 
  ggtitle(label="Rezago Social")

Parecido a los 2 anteriores, nuevamente en el norte entre menor sea el rezago social hay más tendencia a que en los alrededores también sea bajo. Y en este caso, si se pueden notar bien los clusters.

Mapa: LISA Analysis para Densidad Poblacional

lisa_pppe<-local_moran(queen_w, connectivityM["popden2020"]) 

connectivityM1$popden2020<-as.factor(lisa_pppe$GetClusterIndicators())
levels(connectivityM1$popden2020)<-lisa_pppe$GetLabels() 

ggplot(data=connectivityM1) + 
  geom_sf(aes(fill=popden2020)) + 
  ggtitle(label="Densidad de población")

Sin prácticamente clusters o muy pocos de estos, esta variable solamente sigue la tendencia de las 3 gráficas anteriores en pocos municipios.

Mapa: LISA Analysis para Hospitales por Habitantes

lisa_pppe<-local_moran(queen_w, connectivityM["hospitalespHabitantes"]) 

connectivityM1$hospitalespHabitantes<-as.factor(lisa_pppe$GetClusterIndicators())
levels(connectivityM1$hospitalespHabitantes)<-lisa_pppe$GetLabels() 

ggplot(data=connectivityM1) + 
  geom_sf(aes(fill=hospitalespHabitantes)) + 
  ggtitle(label="Hospitales por Habitantes")

Prácticamente sin tendencias, más que en Chiapas donde podemos ver una relación LL, que indica que si hay pocos hospitales para cada 10 mil habitantes en un municipio, los otros también tienden a tener pocos.

Principales hallazgos:

  • Las variables socio-económicas y de acceso a servicios de salud son las mejores para este análisis.
  • Las variables requirieron transformaciones donde: si eran porcentuales se dividieron entre 100 y la de Densidad de Población fue transformada a su valor logarítmico.
  • En las gráficas exploratorias de las variables podemos observar como el sur del país tiende a tener las peores cifras en términos socio-encomios y de acceso a servicios de salud.
  • Y en contra parte, en las gráficas de análisis LISA, el norte suele presentar las autocorrelaciones espaciales positivas LL en variables deseadas, indicando posibles municipios de plusvalía.
  • Por lo tanto, es posible intuir aún sin los modelos de regresión que posiblemente el sur del país fue más afectado por el Covid dada sus condiciones socio-económicas y de acceso a servicios de salud, mientras que el norte probablemente presente mayores contagios y a mayor velocidad.
  • Todas las variables cuentan con autocorrelación espacial positiva.

Conclusiones y preguntas detonadoras: