UN EJEMPLO CON DATOS DEL COVID

Anteriormente, leí la base de datos completa en R. (¿cómo leer una base de datos tan grande? un buen reto). A continuación, realicé un filtrado (Jalisco) y, para fines didácticos y que ustedes puedan reproducir este ejercicio, guardé solo las columnas que me pueden llegar a ser de utilidad.

Pueden no correr las líneas siguientes ya que les envío el archivo datos_jau2.csv

#datos_jau2<-datos_jau %>% dplyr::filter(ENTIDAD_RES=="14")
#write.csv(datos_jau2[,-c(1:4,32:40)],file="datos_jau2.csv")
datos_jau2<-read.csv("datos_jau2.csv")

Cambié el formato de la columna “FECHA_INGRESO”, y obtuve la frecuencia de cada fecha. Ahí agregué una columna con “unos” para que la suma de esa columna sea la frecuencia.

datos_jau2$FECHA_INGRESO<-as.Date(datos_jau2$FECHA_INGRESO,format="%Y-%m-%d")
datos_jau2$dumx<-rep(1,length(datos_jau2$FECHA_INGRESO))

serie_tiempo<-aggregate(dumx~FECHA_INGRESO,data=datos_jau2,sum)
plot(serie_tiempo,type="l")

Un gráfico de Serie de Tiempo interactivo es brindado por la paquetería dygraphs, la cual requiere de ‘input’ una serie de tiempo en formato xts

don <- xts(x=serie_tiempo,order.by=serie_tiempo$FECHA_INGRESO)
dygraph(don) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="#D8AE5A") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)

Ahora, hago el mismo tratamiento para la columna FECHA_DEF, asignando NA a las fechas de los pacientes que no fallecieron. También se genera una gráfica interactiva con Roll-bar.

datos_jau2$FECHA_DEF<-as.Date(datos_jau2$FECHA_DEF,format="%Y-%m-%d")
serie_tiempo2<-aggregate(dumx~FECHA_DEF,FUN = sum,data=datos_jau2)

ser_tiempo1<-data.frame(fechas=seq(serie_tiempo[1,1],serie_tiempo[nrow(serie_tiempo),1],by="day"),casos=rep(0,nrow(serie_tiempo)))

defu_fechas<-merge(ser_tiempo1,serie_tiempo2,by.x="fechas",by.y="FECHA_DEF",sort=TRUE,all.x=TRUE,incomparables = "0")

defu_fechas[which(is.na(defu_fechas$dumx)==TRUE),3]<-0
don2 <- xts(x=defu_fechas$dumx,order.by=defu_fechas$fechas)

dygraph(don2) %>%
  dyOptions(labelsUTC = TRUE, fillGraph=TRUE, fillAlpha=0.1, drawGrid = FALSE, colors="#D8AE5A") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)  %>%
  dyRoller(rollPeriod = 1)

Para comenzar con la graficación, manejaré el ejemplo de las regiones de Jalisco. Ver archivo anexo llamado “Regiones_jalisco.xlsx” Ahora, y a forma de ejercicio, se llamarán los archivos espaciales con otra paquetería sf, la cual, aunque es más limitada en comparación con la utilizada anteriormente, es muy buena para nuestros fines (función st_read). ### cbind

#Regiones Jalisco
regiones_jal<-read_excel("Regiones_jalisco.xlsx")
#Mapa de Jalisco
jal_map<-st_read("jalisco_mun.shp",options = "ENCODING=WINDOWS-1252")
## options:        ENCODING=WINDOWS-1252 
## Reading layer `jalisco_mun' from data source `C:\Users\Hector de la Torre\Documents\CIDE\Progra\finalsR\jalisco_mun.shp' using driver `ESRI Shapefile'
## Simple feature collection with 125 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -105.6954 ymin: 18.92587 xmax: -101.5105 ymax: 22.75025
## CRS:            4326
#Le agrego una columna correspondiente a la región al archivo espacial y le llamo ```jal_map2```
jal_map2<-sp::merge(x=jal_map,y=regiones_jal,by.x="NOMGEO",by.y="Nombre",all.x=TRUE,sort=FALSE)

Grafico de forma estática dicho mapa utilizando la paquetería ggplot2

jal_map2 %>% ggplot() + geom_sf(aes(fill=Region))+labs(fill="Región")

Genero algunos estadísticos útiles a partir de la base de datos:

#Casos positivos Jalisco
casos_mun<-aggregate(dumx~MUNICIPIO_RES,dat=datos_jau2,FUN=sum)

#Fallecimientos
casos_defun<-which(is.na(datos_jau2$FECHA_DEF)!=TRUE)
falle_mun<-aggregate(dumx~MUNICIPIO_RES,dat=datos_jau2[casos_defun,],FUN=sum)

Vuelvo a utilizar los datos de los resultados del INEGI para obtener la población total de cada municipio según el Censo de Población y Vivienda 2020.

#Población jalisco
censo_jal<-read_excel("ITER_14XLSX20.xlsx")
censo_jal_mun<-censo_jal %>% filter(LOC=="0000")

Comienzo a unir al archivo espacial la información sobre población total, casos COVID y fallecidos. Después, obtengo el índice de positivos en la población y porcentajes de fallecimiento.

#Mapa + población
jal_map3<-merge(jal_map2,censo_jal_mun[,c(3,10)],by.x="CVE_MUN",by.y="MUN",all.x=TRUE,sort=FALSE)
#Mapa + casos
library(stringr)
casos_mun$MUNICIPIO_RES<-str_pad(casos_mun$MUNICIPIO_RES, 3, side="left", pad="0")
jal_map4<-merge(jal_map3,casos_mun,by.x="CVE_MUN",by.y="MUNICIPIO_RES",all.x=TRUE,sort=FALSE)
colnames(jal_map4)[7]<-"CASOS_POS"
#Mapa + fallecidos
falle_mun$MUNICIPIO_RES<-str_pad(falle_mun$MUNICIPIO_RES, 3, side="left", pad="0")
jal_map5<-merge(jal_map4,falle_mun,by.x="CVE_MUN",by.y="MUNICIPIO_RES",all.x=TRUE,sort=FALSE)
colnames(jal_map5)[8]<-"FALLE"
#jal_map5$FALLE[125]<-0 #Hay un municipio donde no hubo fallecidos

jal_map5$Pro_positi<-jal_map5$CASOS_POS/as.numeric(jal_map5$POBTOT)
jal_map5$Pro_falle<-jal_map5$FALLE/jal_map5$CASOS_POS

Primero, generaré dos mapas estáticos para estos dos índices Fallecimientos

jal_map5 %>% ggplot() + geom_sf(aes(fill=Pro_falle), size = 0.25)+ 
  scale_fill_distiller(name="% Falle", palette = "YlGn", breaks = pretty_breaks())

Casos positivos

jal_map5 %>% ggplot() + geom_sf(aes(fill=Pro_positi))+ 
  scale_fill_distiller(name="% Positivos", palette = "Reds", breaks = pretty_breaks())

Utilizaré estos dos índices para generar dos mapas vivos. Ojo, aquí agrego ya leyenda de la información y utilizo algo llamado pretty_breaks() de la paquetería scales para genera los breaks de manera automática.

Para el porcentaje de positivos:

cortess <- c(0,0.01532,0.02233,0.03166,Inf)
colores <- colorBin( palette="Reds", domain=jal_map5$Pro_positi, na.color="transparent", bins=cortess)

textoss <- paste(
  "Municipio : ",jal_map5$NOMGEO,"<br/>",
  "% Positivos: ", round((jal_map5$Pro_positi)*100,2),"<br/>", 
  "% Fallecimientos: ", round(jal_map5$Pro_falle*100,2) )  %>% lapply(htmltools::HTML)

  leaflet(data=jal_map5) %>% 
    addTiles() %>% 
    addPolygons(label = textoss,
                fillColor = colores(jal_map5$Pro_positi),
                fillOpacity = 0.9,color="#EFC000FF",weight=2) %>% 
    leaflet::addLegend("bottomright", pal = colores, values = ~Pro_positi,
              title = "% Positivos",
              #labFormat = labelFormat(prefix = "$"),
              opacity = 1 )

Para el porcentaje de fallecimientos:

#% Fallecimientos
library(leaflet)
cortess <- c(0,0.05719,0.08206,0.11085,Inf)
colores <- colorBin( palette="Blues", domain=jal_map5$Pro_falle, na.color="transparent", bins=cortess)

leaflet(data=jal_map5) %>% 
  addTiles() %>% 
  addPolygons(label = textoss,
              fillColor = colores(jal_map5$Pro_falle),
              fillOpacity = 0.9,color="#EFC000FF",weight=2) %>% 
  leaflet::addLegend("bottomright", pal = colores, values = ~Pro_falle,
            title = "% Fallecimientos",
            #labFormat = labelFormat(prefix = "$"),
            opacity = 1 )

Lo anterior está disponible en (https://rpubs.com/coronamexico/842305)