Introducción

El control de la expansión del coronavirus Sars-COV-2 presenta serios desafíos para los sistemas sanitarios de todo el mundo, por lo que resulta importante el poder identificar y cuantificar la afectacción de las distintas zonas o niveles geo-políticos de manera a tomar decisiones adecuadas teniendo en cuenta cada situación en particular (Población y estado del sistema sanitario a nivel local principalmente).

Es así que se plantea el problema del análisis y reporte automatizado con un enfoque geográfico, de manera a contar con información cuantificada en diferentes niveles que reflejen el desarrollo de la pandemia a nivel local (relativamente). En ese sentido, se ha desarrollado este script tratando de proveer la mayor cantidad de información relevante a cada escala de trabajo, en el caso de Paraguay.

Obbjetivos principales:

  • Cuantificar e ilustrar la cantidad de casos confirmados acumulados y la cantidad de casos activos nivel de departamentos y distritos, así también se presentan datos como el total de ambas variables a la fecha de actualización y un gráfico de barras presentando de manera descendente las unidades geográficas con mayor cantidad de casos confirmados

  • Cuantificar y visualizar la variación de movilidad en Paraguay a través de departamentos

Consideraciones:

  • Los datos trabajados son en última instancia correspondientes a los casos comunitarios o fuera de albergues (Derivados del contagio por transmisión comunitaria).
  • Los insumos necesarios son:
    • Datos tabulados de casos de COVID-19
    • Archivos shapefile correspondientes a las unidades geográficas a mapear (Distritos y departamentos).
    • Datos de movilidad de Google
  • Se basa principalmente en tidyverse, utilizando ggplot2 para la composición de los mapas
  • Cada mapa fue elaborado teniendo en cuenta la siguiente configuración de impresión (width = 18, height = 11, dpi = 300).

Metodología

El flujo de trabajo fue le siguiente:

  1. Carga de librerias
  2. Directorio de trabajo
  3. Importación los datos a utilizar
  4. Preparación datos geográficos
  5. Procesamiento de los datos COVID (Renombrar variables, cambiar formatos de campo, asignaciones, categorización, operaciones aritméticas, agrupaciones según conjuntos de variables).
  6. Composición de mapas
  7. Exportación

Carga de librerias a utilizar

##Library

  library(sf) # Manejo de datos vectoriales
  library(leaflet) # Composición de mapas interactivos
  library(RColorBrewer) # Paletas
  library(tidyverse) # Procesamiento datos
  library(zoo) # Rolling functions
  library (tmaptools) ## Palette_explorer()
  library(readxl) ## Importar datos excel
  library(htmlwidgets) # Guardado de archivos .html
  library(ggplot2) # Visualizaciones, mapas
  library(ggspatial) # No
  library(ggrepel)
  library(scales)
  library(plotly)
  library(ggcorrplot)

  # Covid
  library(covid19mobility)
  library(COVID19)

Establecimiento del directorio de trabajo

Importación de datos

Aquí se importan los datos necesarios, es decir los datos vectoriales en formato “shapefile” y las planillas “csv” correspondientes a los datos de COVID-19 en Paraguay. En relación a los datos shapefile es importante mencionar que, estos fueron previamente procesados, corrigiéndose la topología y la tabla de atributos asociada (Códigos de departamentos y distritos), los datos correspondientes a departamentos y distritos fueron de esta manera adaptados de (DGEEC,2012). En cuanto a las zonas epidemiológicas estas fueron generadas agregando los departamentos que cada una abarca. De la misma manera, en orden a generar mapas de burbujas, se extrajeron los centroides de cada unidad geográfica y a estos le fue asignada toda la información espacial correspondiente.
# Se creo un objeto espacial el cual corresponde a los datos vectoriales de departamentos
dpto_shape <-
  read_sf('SHAPEFILES/wgs84_projected/covid_py_dptos.shp',
          stringsAsFactors = FALSE)
# Se creo un objeto espacial el cual corresponde a los datos vectoriales de distritos
# Además, se agrega un campo conteniendo el código del Dpto. y distrito de manera a identificar cada distrito de manera única
dist_shape <-
  read_sf('SHAPEFILES/wgs84_projected/covid_py_dist.shp',
          stringsAsFactors = FALSE) %>% mutate(dpto_dist_cod = paste(dpto_cod, dist_cod, sep = '_'))
##Capas de puntos
# Objeto espacial de puntos correspondiente a los centroides de los polígonos del archivo shapefile de departamentos
dpto_points <-
  read_sf('SHAPEFILES/wgs84_projected/covid_py_dpto_points.shp',
          stringsAsFactors = FALSE)%>% mutate(dpto_desc = ifelse(dpto_desc=="ÑEEMBUCU", "ÑEEMBUCU", dpto_desc))
# Objeto espacial de puntos correspondiente a los centroides de los polígonos del archivo shapefile de distritos
# Además, se agrega un campo conteniendo el código del Dpto. y distrito de manera a identificar cada distrito de manera única
dist_points <-
  read_sf('SHAPEFILES/wgs84_projected/covid_py_dist_points.shp',
          stringsAsFactors = FALSE) %>% mutate(dpto_dist_cod = paste(dpto_cod, dist_cod, sep = '_'))
# Objeto conteniendo los registros de casos de COVID-19 en Paraguay obtenidos del Ministerio de Salud Pública y Bienestar Social
# Aquí siempre es importante verificar el archivo de entrada, pues el mismo va siendo modificado en términos de formato por la misma fuente (MSyBS), principalemte separadores y formato de fecha
datos_brutos <-
  data.frame(
    read.csv(
      "Descargar_datos_Datos_completos_data.csv",
      sep = ';',
      encoding = "UTF-8",
      stringsAsFactors = FALSE
    )
  )
# Estos son datos globales de movilidad reportados por Google, los mismos se descargan de https://www.google.com/covid19/mobility/
global_movility <-
  data.frame(read.csv("Global_Mobility_Report.csv", sep = ',', encoding = "UTF-8", stringsAsFactors = FALSE))

#A partir de aquí los datos son directamente importados de un repositorio, dado que en los objetos anteriores se cargaban los datos previamente descargados de manera local

# Estos son datos de covid19 descargados de un repositorio, los mismos contienen 36 variables, de estas se usa la cantidad de casos acumulados, muertes y recuperados ya que en los datos
# descargados de la página de reporte del MSyBS no se cuenta con esas variebles según registro, hay que mencionar que estos datos ya se encuentran trabajados
Paraguay <- covid19(
  country = "PY",
  level = 1,
  start = "2010-01-01",
  end = Sys.Date(),
  raw = TRUE,
  vintage = FALSE,
  verbose = TRUE,
  cache = TRUE,
  wb = NULL,
  gmr = NULL,
  amr = NULL
) %>% mutate(deaths = replace_na(deaths,0), recovered = replace_na(recovered,0)) %>%  mutate(active_c = confirmed - (recovered+deaths))# Se calcula la cantidad de datos para este set de
## We have invested a lot of time and effort in creating COVID-19 Data Hub, please cite the following when using it:
## 
##   Guidotti, E., Ardia, D., (2020), "COVID-19 Data Hub", Journal of Open
##   Source Software 5(51):2376, doi: 10.21105/joss.02376.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Article{,
##     title = {COVID-19 Data Hub},
##     year = {2020},
##     doi = {10.21105/joss.02376},
##     author = {Emanuele Guidotti and David Ardia},
##     journal = {Journal of Open Source Software},
##     volume = {5},
##     number = {51},
##     pages = {2376},
##   }
## 
## To retrieve citation and metadata of the data sources see ?covid19cite. To hide this message use 'verbose = FALSE'.
# Descarga los datos de movilidad por fecha de un repositorio de datos mediante el paquete covid19mobility
goog <- refresh_covid19mobility_google_subregions()
## Warning: 15010 parsing failures.
##  row        col           expected                  actual                                                                  file
## 2046 metro_area 1/0/T/F/TRUE/FALSE Kabul Metropolitan Area 'https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv'
## 2047 metro_area 1/0/T/F/TRUE/FALSE Kabul Metropolitan Area 'https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv'
## 2048 metro_area 1/0/T/F/TRUE/FALSE Kabul Metropolitan Area 'https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv'
## 2049 metro_area 1/0/T/F/TRUE/FALSE Kabul Metropolitan Area 'https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv'
## 2050 metro_area 1/0/T/F/TRUE/FALSE Kabul Metropolitan Area 'https://www.gstatic.com/covid19/mobility/Global_Mobility_Report.csv'
## .... .......... .................. ....................... .....................................................................
## See problems(...) for more details.
## Joining, by = c("location_code", "sub_region_1")

Preparación de datos geográficos para cada nivel o escala a mapear

Se procesaron los datos importados de los shapefile, preparándolos para agregarlos luego con los datos de las variables de COVID-19

# Creación de un dataframe alojando los datos de coodernadas de cada unidad a mapear (Los puntos son los centroides de los polígonos de cada distrito, departamento y zona)
 dpto_cords <-
  rename(
    dpto_points,
    dpto_desc = dpto_desc,
    dpto_cod = dpto_cod,
    "lon_dpto" = lon,
    "lat_dpto" = lat
  ) %>% select(dpto_desc, dpto_cod, lon_dpto, lat_dpto)## Se extraen las variables de dpto y dist desc y cod
dist_cords <-
  rename(
    dist_points,
    dist_desc = dist_desc,
    dpto_dist_cod = dpto_dist_cod,
    "lon_dist" = lon,
    "lat_dist" = lat
  ) %>% select(dist_desc, dpto_dist_cod, lon_dist, lat_dist)

Procesamiento de datos con Dplyr

En esta sección se realizan tareas como: Renombrar campos (Variables), transformar los datos de las variables de manera a que sea más sencillo obtener información de ellas a partir de una escala dada.
#Filtro de los datos de movilidad para Paraguay y se crean nuevas variables para poder agregar los datos según fechas y departamentos

# ggplot(Paraguay_mob, aes(fecha, retail_and_recreation_percent_change_from_baseline)) + geom_line(size = 2, color = "blue") visualizacion de caracter
#Renombrando los campos de la base de datos descargada de Salud
datos_brutos <- rename(
  datos_brutos,
  fecha = Fecha.Confirmacion, # Debido a los frecuentes inconvenientes con este campo, es altamente recomendable trabajarlo de manera a eliminar previamente las fuentes posibles de error
  sexo = Sexo,
  dpto_desc = Departamento.Residencia,
  dist_desc = Distrito.Residencia,
  albergue = En.Albergue.,
  dist_cod = Cod.Distrito,
  dpto_cod = Cod.Dpto,
  edad = Edad
)


# Se seleccionan y transforman los datos de manera a poder procesarlos por cada unidad espacial (Distritos, departamentos y zonas epidemiológicas)
# Siempre verificar los datos de fecha principalmente...
datos_covid <- select(datos_brutos,
                      fecha,
                      sexo,
                      dpto_desc,
                      dist_desc,
                      albergue,
                      dist_cod,
                      dpto_cod,
                      edad) %>% mutate(
# Se cammbia el formato del campo fecha, de cadena a Date y se establece que es que (mes, día, año) en el formato específico
                        fecha = as.Date(fecha, "%d/%m/%Y"), # Importante: Verificar el formato de entrada, en caso de que se haya modificado simplemente cambiar la interpretacion segun sea %d: dia, %m: mes, %Y: año
# Se crea un campo para identificar los distritos de cada registro
                        dpto_dist_cod = paste(dpto_cod, dist_cod, sep = '_'),
# Se asigna NA (Not available)
                        edad = ifelse(edad < 120, edad, NA),
                        dpto_desc = ifelse(dpto_desc=="ÑEEMBUCU", "ÑEEMBUCU", dpto_desc),
# Se asigna la zona epidemiológica a la que pertenece cada registro a través de ifelse sucesivos y los códigos de cada Dpto.
                        zona_epi = ifelse(
                          dpto_cod %in% c(0, 11),
                          "Asunción y Central",
                          ifelse(
                            dpto_cod %in% c(3, 4, 5, 6, 10),
                            "Centro Este",
                            ifelse(
                              dpto_cod %in% c(15, 16, 17),
                              "Chaco",
                              ifelse(dpto_cod %in% c(1, 2, 13, 14), "Norte", "Sur")
                            )
                          )
                        )
                      ) %>% mutate(
# Se crea el campo iso_active y se popula como si fuera un indice
                        iso_active = 1:length(fecha),
# Se calcula la semana ISO de confirmación del caso para cada registro
                        iso_week = as.integer(fecha %>% strftime(format = "%V")),
# Si la fecha de confirmación + 21 días (3 semanas) es menor a la fecha de hoy(Fecha de ejecución del script), entonces se asignará un 0 (inactivo), si no asignará un 1 indicando que el registro sigue activo  
                        activo = ifelse((fecha + 21) < Sys.Date(), 0, 1),
# Se calcula la fecha del último día actividad, y se asigna a menos que dicha fecha sea mayor a la fecha del sistema, en cuyo caso se asigna la fecha del sistema
                        fecha_activo = ifelse((fecha + 21) < Sys.Date(),
                                              as.character((fecha + 21)),
                                              as.character(Sys.Date())
                        )
                      ) %>% mutate(
# Se asigna la semana iso de la última fecha de actividad del registro
                        iso_last =  as.integer(strftime(as.Date(fecha_activo), format = "%V")) 
                      )
ind_fases <-
  c(
    "Pre-cuarentena",
    "Cuarentena parcial",
    "Cuarentena total",
    "Flexibilización (de preparación)",
    "Flexibilización (leve)",
    "Flexibilización (Moderada)",
    "Fase 4"
  )

datos_covid_py <-
  datos_covid %>%  filter(datos_covid$fecha <= as.Date(max(Paraguay$date), "%Y-%m-%d") & datos_covid$albergue == "NO") %>% mutate(
    iso_rank = 1 + (iso_last - iso_week),
    fase = ifelse(
      fecha >= "2020-03-11" &
        fecha <= "2020-03-19",
      "Cuarentena parcial",
      ifelse(
        fecha >= "2020-03-20" &
          fecha <= "2020-05-03",
        "Cuarentena total",
        ifelse(
          fecha >= "2020-05-04" &
            fecha <= "2020-05-24",
          "Flexibilización (de preparación)",
          ifelse(
            fecha >= "2020-05-25" &
              fecha <= "2020-06-14",
            "Flexibilización (leve)",
            ifelse(
              fecha >= "2020-06-15" &
                fecha <= "2020-07-19",
              "Flexibilización (Moderada)",
              ifelse(fecha >= "2020-07-20",
                     "Fase 4", "Pre-cuarentena")
            )
          )
        )
      )
    )
  ) 

# Se agregan los registros de cada caso por fecha y se cuantifica la cantidad de casos por fecha y el acumulado a cada fecha
datos_por_dia <-datos_covid_py %>%  group_by(fecha) %>% summarise(conf_dia = n()) %>% mutate(acum = cumsum(conf_dia))

# Los registros de covid se agrupan por fecha y departamento y se cuenta la cantidad de confirmados y activos
byDateDpto_covidpy_df <- datos_covid_py %>% group_by(fecha, dpto_desc) %>% summarise("cant_casos" = n(), "activo" = sum(as_tibble(activo))) %>% mutate(identificador = paste(fecha,dpto_desc, sep = "_"))
datos_e1 <- #Datos actuales, es decir un set de datos sin acumular, ser utilizados para ilustrar la situacion actual
  mutate(
    datos_covid_py,
    #Se crea la variable 'hombres' en donde a través de una función ifelse se evalua si el valor de los registros de la variable 'sexo' es MASCULINO o no, en caso de que sí lo sea se asigna un 1 y en caso contrario un 0
    hombres = ifelse(datos_covid_py$sexo == "MASCULINO", 1, 0),
    #Se crea la variable 'mujeres' en donde a través de una función ifelse se evalua si el valor de los registros de la variable 'sexo' es FEMENINO o no, en caso de que sí lo sea se asigna un 1 y en caso contrario un 0
    mujeres = ifelse(datos_covid_py$sexo == "FEMENINO", 1, 0),
    #Se crea una variable 'menor_40' en donde a través de una función ifelse se evalua si el valor de los registros de la variable 'edad' es menor a 40, en caso de que sí lo sea se asigna un 1 y en caso contrario (Si es igual o mayor) un 0
    menor_40 = ifelse(datos_covid_py$edad < 40, 1, 0),
    #Se crea una variable 'entre_40_50' en donde a través de una función ifelse se evalua si el valor de los registros de la variable 'edad' es mayor o igual a 40 y menor o igual a 50, en caso de cumplir dichas condiciones se asigna un 1 y en caso contrario se asigna 0
    entre_40_50 = ifelse(datos_covid_py$edad >= 40 &
                           datos_covid_py$edad <= 50, 1, 0),
    #Se crea una variable 'mayor_50' en donde a través de una función ifelse se evalua si el valor de los registros de la variable 'edad' es mayor a 50, en caso de cumplir dicha condicion se asigna un 1 y en caso contrario se asigna 0
    mayor_50 = ifelse(datos_covid_py$edad > 50, 1, 0),
    #Se crea el campo dpto_dist_cod a fin de asignar a cada registro un código constituido por el cod del departamento mas el codigo de distrito
    dpto_dist_cod = paste(dpto_cod, dist_cod, sep = '_')
  )
#Transformacion de variables con valores 'cadena' a numericos, de manera a resumirlos luego a través de operaciones aritméticas
#Resumen de variables por unidad geográfica (Datos actuales)
datos_dptos <- # Resumen por departamentos----
datos_e1 %>% group_by(dpto_cod) %>% summarise(
  ##Se agrupan los registros por departamentos en este caso 12 departamentos
  "cant_casos" = n(),
  ##Se cuenta el número de registros de agrupación
  "activo" = sum(as_tibble(activo)),
  ## Suma de de esta columna para cada agrupación
  "mujeres" = sum(as_tibble(mujeres)),
  "hombres" = sum(as_tibble(hombres)),
  "menor_40" = sum(as_tibble(menor_40), na.rm = TRUE),
  "entre_40_50" = sum(as_tibble(entre_40_50), na.rm = TRUE),
  "mayor_50" = sum(as_tibble(mayor_50), na.rm = TRUE)
) %>% mutate(
  pred_sex = ifelse(
    ## Se añade un campo "pred_sex" y se calcula cual es el sexo "más afectado" en esa
    hombres < mujeres,
    ## localidad, teniendo en cuenta los campos calculados anteriormente
    "Femenino",
    ifelse(mujeres == hombres, "Ambos", "Masculino")
  ),
  cant_casos = replace_na(cant_casos, 0),
  etario_pred = ifelse(
    ## Se añade un campo "etario_pred" y se calcula cual es el sexo "más afectado" en esa
    menor_40 > entre_40_50 & ## unidad geográfica
      menor_40 > mayor_50,
    "Menor a 40",
    ifelse(
      menor_40 == entre_40_50 &
        menor_40 > mayor_50,
      "Menor a 50",
      ifelse(
        entre_40_50 > mayor_50 &
          entre_40_50 > menor_40,
        "Entre 40 y 50",
        ifelse(
          entre_40_50 == mayor_50,
          "Mayor a 40",
          ifelse(mayor_50 == menor_40, "Menor a 40 y mayor a 50", "Mayor a 50")
        )
      )
    )
  )
) %>% merge(dpto_cords, by = 'dpto_cod') %>%  select(-geometry)## Se agregan las coordenadas de centroides de cada departamento

## Resumen por distritos----
datos_dist <-
  datos_e1 %>% group_by(dpto_dist_cod) %>% summarise(
    "cant_casos" = n(),
    "activo" = sum(as_tibble(activo)),
    "mujeres" = sum(as_tibble(mujeres)),
    "hombres" = sum(as_tibble(hombres)),
    "menor_40" = sum(as_tibble(menor_40), na.rm = TRUE),
    "entre_40_50" = sum(as_tibble(entre_40_50), na.rm = TRUE),
    "mayor_50" = sum(as_tibble(mayor_50), na.rm = TRUE)
  ) %>% mutate(
    pred_sex = ifelse(
      hombres < mujeres,
      "Femenino",
      ifelse(mujeres == hombres, "Ambos", "Masculino")
    ),
    cant_casos = replace_na(cant_casos, 0),
    etario_pred = ifelse(
      menor_40 > entre_40_50 &
        menor_40 > mayor_50,
      "Menor a 40",
      ifelse(
        menor_40 == entre_40_50 &
          menor_40 > mayor_50,
        "Menor a 50",
        ifelse(
          entre_40_50 > mayor_50 &
            entre_40_50 > menor_40,
          "Entre 40 y 50",
          ifelse(
            entre_40_50 == mayor_50,
            "Mayor a 40",
            ifelse(mayor_50 == menor_40, "Menor a 40 y mayor a 50", "Mayor a 50")
          )
        )
      )
    )
  ) %>% merge(dist_cords, by = 'dpto_dist_cod') %>% select(-geometry)

#Resumen de datos por unidad geográfica y segun fases (Datos actuales por fases)
##Datos fases departamentos----
datos_dpto_fases <-
  datos_e1 %>% group_by(dpto_cod, fase) %>% summarise(
    ##Se agrupan los registros por departamentos en este caso 12 departamentos
    "cant_casos" = n(),
    ##Se cuenta el número de registros de agrupación
    "activo" = sum(as_tibble(activo)),
    ## Suma de de esta columna para cada agrupación
    "mujeres" = sum(as_tibble(mujeres)),
    "hombres" = sum(as_tibble(hombres)),
    "menor_40" = sum(as_tibble(menor_40), na.rm = TRUE),
    "entre_40_50" = sum(as_tibble(entre_40_50), na.rm = TRUE),
    "mayor_50" = sum(as_tibble(mayor_50), na.rm = TRUE)
  ) %>% mutate(
    pred_sex = ifelse(
      ## Se añade un campo "pred_sex" y se calcula cual es el sexo "más afectado" en esa
      hombres < mujeres,
      ## localidad, teniendo en cuenta los campos calculados anteriormente
      "Femenino",
      ifelse(mujeres == hombres, "Ambos", "Masculino")
    ),
    cant_casos = replace_na(cant_casos, 0),
    etario_pred = ifelse(
      ## Se añade un campo "etario_pred" y se calcula cual es el sexo "más afectado" en esa
      menor_40 > entre_40_50 & ## unidad geográfica
        menor_40 > mayor_50,
      "Menor a 40",
      ifelse(
        menor_40 == entre_40_50 &
          menor_40 > mayor_50,
        "Menor a 50",
        ifelse(
          entre_40_50 > mayor_50 &
            entre_40_50 > menor_40,
          "Entre 40 y 50",
          ifelse(
            entre_40_50 == mayor_50,
            "Mayor a 40",
            ifelse(mayor_50 == menor_40, "Menor a 40 y mayor a 50", "Mayor a 50")
          )
        )
      )
    )
  ) %>% merge(dpto_cords, by = 'dpto_cod') %>%  select(-geometry) %>% mutate(fase = factor(fase, levels = ind_fases))

##Datos fases distritos----
datos_dist_fases <-
  datos_e1 %>% group_by(dpto_dist_cod, fase) %>% summarise(
    ##Se agrupan los registros por departamentos en este caso 12 departamentos
    "cant_casos" = n(),
    ##Se cuenta el número de registros de agrupación
    "activo" = sum(as_tibble(activo)),
    ## Suma de de esta columna para cada agrupación
    "mujeres" = sum(as_tibble(mujeres)),
    "hombres" = sum(as_tibble(hombres)),
    "menor_40" = sum(as_tibble(menor_40), na.rm = TRUE),
    "entre_40_50" = sum(as_tibble(entre_40_50), na.rm = TRUE),
    "mayor_50" = sum(as_tibble(mayor_50), na.rm = TRUE)
  ) %>% mutate(
    pred_sex = ifelse(
      ## Se añade un campo "pred_sex" y se calcula cual es el sexo "más afectado" en esa
      hombres < mujeres,
      ## localidad, teniendo en cuenta los campos calculados anteriormente
      "Femenino",
      ifelse(mujeres == hombres, "Ambos", "Masculino")
    ),
    cant_casos = replace_na(cant_casos, 0),
    etario_pred = ifelse(
      ## Se añade un campo "etario_pred" y se calcula cual es el sexo "más afectado" en esa
      menor_40 > entre_40_50 & ## unidad geográfica
        menor_40 > mayor_50,
      "Menor a 40",
      ifelse(
        menor_40 == entre_40_50 &
          menor_40 > mayor_50,
        "Menor a 50",
        ifelse(
          entre_40_50 > mayor_50 &
            entre_40_50 > menor_40,
          "Entre 40 y 50",
          ifelse(
            entre_40_50 == mayor_50,
            "Mayor a 40",
            ifelse(mayor_50 == menor_40, "Menor a 40 y mayor a 50", "Mayor a 50")
          )
        )
      )
    )
  ) %>% merge(dist_cords, by = 'dpto_dist_cod') %>%  select(-geometry) %>% mutate(fase = factor(fase, levels = ind_fases))

Resumen de variables por unidad geográfica y unión

En esta sección los datos transformados se agrupan y cuantifican por unidad geográfica (Departamentos y distritos), así también se determina la predominancia de ciertas variables en cada localidad. Las variables trabajadas son:

  • Cantidad de casos
  • Cantidad de casos activos

Finalmente a cada objeto creado se agregan los campos correspondientes a las coordenadas de los centroides de polígonos de las unidades geográficas

## Se realiza un join entre los sets procesados y los archivos espaciales----
  ##Poligonos
  dpto_shape1 <-
    merge(dpto_shape, datos_dptos, by = 'dpto_cod', all.x = TRUE) %>% mutate(cant_casos = replace_na(cant_casos, 0)) %>% rename(
      dpto_desc = dpto_desc.x)
  
  dist_shape1 <-
    merge(dist_shape, datos_dist, by = "dpto_dist_cod", all.x = TRUE) %>% rename(dist_desc = dist_desc.x)
  
  ##Puntos
  dpto_points1 <-
    merge(dpto_points, datos_dptos, by = 'dpto_cod', all.x = TRUE)
  dist_points1 <-
    merge(dist_points, datos_dist, by = "dpto_dist_cod", all.x = TRUE)

Procesamiento de datos de movilidad

Se filtran los datos de movilidad para Paraguay, luego se asigna un código de departamento según DGEEC. 20112, y se establece un campo de fecha

## Warning: 229 parsing failures.
## row col expected actual
## 685  -- a number      -
## 686  -- a number      -
## 687  -- a number      -
## 688  -- a number      -
## 689  -- a number      -
## ... ... ........ ......
## See problems(...) for more details.
## Warning in write.csv2(Paraguay_mob, "Paraguay_mob.csv", sep = ";"): attempt to
## set 'sep' ignored
## Warning: 1398 parsing failures.
##  row col expected actual
## 2779  -- a number      -
## 2780  -- a number      -
## 2781  -- a number      -
## 2782  -- a number      -
## 2783  -- a number      -
## .... ... ........ ......
## See problems(...) for more details.

Composición de mapas con R

En esta etapa, los datos preparados y procesados anteriormente serán visualizados a través de mapas estáticos desarrollados a través del paquete Ggplot2.

Composición de mapas base

#ELEMENTOS COMPOSICIÓN DE MAPAS
  ##Mapas base poligonos
  
mapa_base_dpto <- dpto_shape %>%
  ggplot() + geom_sf(fill = "gray90",
                     alpha = 0.8,
                     colour = "Black")

mapa_base_dist <- dist_shape %>%
  ggplot() + geom_sf(fill = "gray90",
                     alpha = 0.1,
                     colour = "Black")

##Mapas base poligonos para burbujas

mapa_base_dpto_burbujas <- dpto_shape %>%
  ggplot() + geom_sf(fill = "gray75",
                     alpha = 0.8,
                     colour = "gray90")
mapa_base_dist_burbujas <- dist_shape %>%
  ggplot() + geom_sf(fill = "gray75",
                     alpha = 0.8,
                     colour = "gray90")

Breaks y labels

#Breaks
activo_breaks <- c(1, 50, 500, 1000, 5000, 10000) # Breaks para la variable "activos"
cant_casos_breaks <-  c(50, 500, 1000, 5000, 10000) 
mov_breaks <- c(-100, -80,-60,-40,-20,0,20,40,60,80,100)
#Labels
mov_labels <- c("-100%", "-80","-60%","-40%","-20%","0%","20%","40%","60%","80%","100%")
cant_casos_labels <- c("1 - 49", "50 - 499", "500 - 999", "1000 - 4999", "5000 - 9999","10000+")
activo_labels <-
  c("1 - 49", "50 - 499", "500 - 999", "1000 - 4999", "5000 - 9999","10000+")# Labels para la variable "activos"
#Fuentes
  titulo_dist <-  "Distribución distrital de casos de COVID-19 en Paraguay"
  titullo_dpto  <- "Distribución departamental de casos de COVID-19 en Paraguay"
  subtitulo <-  paste("Total confirmados:",format(sum(datos_dist$cant_casos, na.rm = TRUE), big.mark = ".", decimal.mark = ","), " y Total activos:",format(sum(datos_dist$activo, na.rm = TRUE), big.mark = ".", decimal.mark = ","), " fuera de albergue")
  fuentes_facet <- 'Fuentes: Dirección General de Vigilancia de la Salud (D.G.V.S.) \n Obs: El análisis se realizó en base a casos fuera de albergue, considerandolos como activos durante 21 días (3 semanas)'
  fuentes_facet_activos <- 'Fuentes: Dirección General de Vigilancia de la Salud (D.G.V.S.). Obs: El análisis se realizó en base a casos fuera de albergue,considerandolos como activos durante 21 días (3 semanas)'
  fuentes <- paste('Fuentes: Datos de https://www.mspbs.gov.py/reporte-covid19.html',"| Datos actualizados al:",strftime(max(datos_covid_py$fecha),format = "%d/%m/%Y"),"\nContacto: Ing. Carlos Giménez: charlieswall@gmail.com / Dr. Pastor Pérez: peperez.estigarribia@pol.una.py")
  fuentes_cant_casos <- 'Fuentes: Dirección General de Vigilancia de la Salud (D.G.V.S.) \n Obs: El análisis se realizó en base a casos fuera de albergue'
  fuentes_sex_pred <- 'Fuentes: Dirección General de Vigilancia de la Salud (D.G.V.S.) \n Obs: El análisis se realizó en base a casos fuera de albergue. Las etiquetas corresponden a la cantidad de confirmados por unidad geográfica'
  fuentes_etario_pred <- 'Fuentes: Dirección General de Vigilancia de la Salud (D.G.V.S.) \n Obs: El análisis se realizó en base a casos fuera de albergue. Las etiquetas corresponden a la cantidad de confirmados por unidad geográfica\n La segunda columna de rangos corresponde a aquellas localidades en donde 2 rangos tienen la misma cantidad de casos'
  
  fuentes_mov <-
  paste(
    'Fuentes: Datos de https://www.google.com/covid19/mobility/ |',
    "Modelo aditivo generalizado (GAM)",
    "\nContacto: Ing. Carlos Giménez: charlieswall@gmail.com / Dr. Pastor Pérez: peperez.estigarribia@pol.una.py"
  )
title_mov <- "Variación de la movilidad en Paraguay debido al COVID-19"
subtitle_mov <- paste("Datos actualizados al:",
                      strftime(max(goog_py$fecha), format = "%d/%m/%Y"))

fases_fechas <- c(
  as.Date("2020-03-11", format = "%Y-%m-%d"),
  as.Date("2020-03-20", format = "%Y-%m-%d"),
  as.Date("2020-05-04", format = "%Y-%m-%d"),
  as.Date("2020-05-25", format = "%Y-%m-%d"),
  as.Date("2020-06-15", format = "%Y-%m-%d"),
  as.Date("2020-07-20", format = "%Y-%m-%d")
)
fases_nombres <- c(
  "Cuarentena Parcial",
  "Cuarentena Total",
  "Flexibilización (de preparación)",
  "Flexibilización (leve)",
  "Flexibilización (moderada)",
  "Flexibilización (moderada x zona)"
)
fases_mov <- data.frame(fases_fechas,fases_nombres) 

###Gráficos de una página

Se elaboran gráficos a ser presentados de manera individual

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 row(s) containing missing values (geom_path).

Variación del desplazamiento a lugares y permanencia en hogares a través del tiempo en paraguay

Gráfico de líneas mostrando la variación del desplazamiento a lugares y permanencia en hogares a través del tiempo en paraguay

plot_py_mov_unificado <-# Gráfico de líneas mostrando la variación del desplazamiento a lugares y permanencia en hogares a través del tiempo en paraguay----
  ggplot(goog_py, aes(x = fecha, y = value)) + geom_line(alpha = 0.15, color = "gray70") + geom_smooth(se = FALSE,
                                                                                                       aes(
                                                                                                         x = fecha,
                                                                                                         y = value,
                                                                                                         color = lugar
                                                                                                       ))  +
  scale_y_continuous(
    name = "Variación porcentual de movilidad",
    breaks = c(-100,-80,-60,-40,-20, 0, 20, 40, 60, 80, 100),
    limits = c(-100, 100)
  ) + scale_x_date(name = "Fecha",
                   date_breaks = "1 month",
                   date_labels = "%b %d") + labs(
                     color = "Tendencia de movilidad",
                     caption = fuentes_mov,
                     title = title_mov,
                     subtitle = subtitle_mov
                   ) +
  theme(
    panel.background = element_blank(),
    panel.grid.major.y = element_line(colour = "Gray75"),
    panel.grid.major.x = element_line(colour = "Gray75"),
    plot.title = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 20,
      face = "bold"
    ),
    plot.subtitle = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 15,
      face = "bold"
    ),
    plot.caption = element_text(
      hjust = 0,
      size = 12,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.position = "right",
    legend.justification = c(0, 1),
    legend.text = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 14
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    axis.title.x = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 14,
      face = "bold"
    ),
    axis.title.y = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 14,
      face = "bold"
    )
  ) +
  annotate(geom = "text",
      x = fases_mov$fases_fechas,
      y = -99,
      label = fases_mov$fases_nombres,
      colour = "darkgray",
      angle = 90,
      hjust = 0
  )  +
annotation_custom(
    plot_curva_acu_acti,
    xmin = 20,
    xmax = -100,
    ymin = -24,
    ymax = -22
  )


 
plot_py_mov_unificado
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 5377 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 row(s) containing missing values (geom_path).

Gráfico de variación de movilidad por cateogría individual

 # Faceted 

plot_py_mov_individual <- ggplot(goog_py,
       aes(x = fecha, y = value)) + geom_line(alpha = 0.15, color = "gray65") + scale_y_continuous(
         breaks = c(-100,-80,-60,-40,-20, 0, 20, 40, 60, 80, 100),
         limits = c(-100, 100)
       ) + scale_x_date(date_breaks = "1 month", date_labels = "%b") +
  facet_wrap(~ lugar)  +
  geom_smooth(se = FALSE,
              aes(x = fecha,
                  y = value,
                  color = lugar)) + geom_hline(yintercept = 0, colour = "red") + labs(color = "Tendencia de movilidad", x = "Fecha", y = "Variación porcentual de movilidad")+ 
   theme(
    panel.background = element_blank(),
    panel.grid.major.y = element_line(colour = "Gray75"),
    panel.grid.major.x = element_line(colour = "Gray75"),
    plot.title = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 20,
      face = "bold"
    ),
    plot.subtitle = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 15,
      face = "bold"
    ),
    plot.caption = element_text(
      hjust = 0,
      size = 12,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.position = "right",
    legend.justification = c(0, 1),
    legend.text = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 14
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    axis.title.x = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 14,
      face = "bold"
    ),
    axis.title.y = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 14,
      face = "bold"
    )
  ) +
annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  )
plot_py_mov_individual
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 5377 rows containing non-finite values (stat_smooth).

Gráficos a insertar dentro de los mapas

Se elaboran los gráficos ques serán insertados en cada mapa

Mapas estáticos con Ggplot2

Mapa de total de casos confirmados por departamento

#Mapa de cantidad de casos confirmados acumulados sin activos por polígonos
es_cant_confirmados_dpto <-
mapa_base_dpto + geom_sf(
  data = dpto_shape1,
  aes(fill = cant_casos),
  colour = "grey75"
) +
  scale_fill_fermenter(
    show.limits = TRUE,
    breaks = cant_casos_breaks,
    na.value = "grey60",
    palette = "YlOrRd",
    direction = 1,
    limits = c(1, max(dpto_shape1$cant_casos, na.rm = TRUE))
  ) +
  geom_label_repel(
    data = datos_dptos,
    aes(x = lon_dpto, y = lat_dpto, label = cant_casos),
    fontface = "bold",
    size = 4,
    alpha = 0.5,
    segment.color = "black",
    force = 0,
    nudge_x = ifelse(
      datos_dptos$dpto_cod == "0",-0.4,
      ifelse(datos_dptos$dpto_cod == "11",-0.1, 0)
    ),
    nudge_y = 0
  ) +
  annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  )+
  annotate("rect", xmin = -63, xmax =-59.7, ymin = -26.1, ymax = -25.3, alpha = .2) +
  annotate("Text", x=-62.9, y=-25.7, label=  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)", size=5,hjust = 0, family = "serif", color = "#4e4d47") +
 annotation_custom(
    bar_cant_casos_dpto,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  ## Para usar los brewer colors se implementa fill_fermenter o fill_distiller
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(fill = "Total confirmados",
       caption = fuentes,
       title = titullo_dpto,
       subtitle =  subtitulo)  +
  guides(
    fill = guide_bins(
      title.position = 'top',
      keywidth = unit(1.5, units = "cm"),
      keyheight = unit(0.8, units = "cm"),
      axis.linewidth = unit(0.4, units = "cm"),
      order = 1,
      show.limits = TRUE
    )
  ) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    panel.grid.major = element_line(size = 0.5, linetype = "solid", colour = "gray80"),
    plot.title = element_text(family = "serif",color = "#4e4d47", size = 20,  face = "bold"),
    plot.subtitle = element_text(family = "serif",color = "#4e4d47", size = 18),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = "right",
    legend.justification = c(0, 1),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "vertical",
  )

es_cant_confirmados_dpto

Mapa de casos confirmados por distrito

# Gráfico de barras para el mapa de casos confirmados por distrito
#Mapa de cantidad de casos confirmados acumulados sin activos por polígonos
es_cant_confirmados_dist <-
mapa_base_dpto + geom_sf(
  data = dist_shape1,
  aes(fill = cant_casos),
  colour = "grey75"
) +
  scale_fill_fermenter(
    show.limits = TRUE,
    breaks = cant_casos_breaks,
    na.value = "grey60",
    palette = "YlOrRd",
    direction = 1,
    limits = c(1, max(dist_shape1$cant_casos, na.rm = TRUE))
  ) +
  annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  )+
  annotate("rect", xmin = -63, xmax =-59.7, ymin = -26.1, ymax = -25.3, alpha = .2) +
  annotate("Text", x=-62.9, y=-25.7, label=  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)", size=5,hjust = 0, family = "serif", color = "#4e4d47") +
 annotation_custom(
    bar_cant_casos_dist,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  ## Para usar los brewer colors se implementa fill_fermenter o fill_distiller
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(fill = "Total confirmados",
       caption = fuentes,
       title = titulo_dist,
       subtitle =  subtitulo)  +
  guides(
    fill = guide_bins(
      title.position = 'top',
      keywidth = unit(1.5, units = "cm"),
      keyheight = unit(0.8, units = "cm"),
      axis.linewidth = unit(0.4, units = "cm"),
      order = 1,
      show.limits = TRUE
    )
  ) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    panel.grid.major = element_line(size = 0.5, linetype = "solid", colour = "gray80"),
    plot.title = element_text(family = "serif",color = "#4e4d47", size = 20,  face = "bold"),
    plot.subtitle = element_text(family = "serif",color = "#4e4d47", size = 18),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = "right",
    legend.justification = c(0, 1),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "vertical",
  )

es_cant_confirmados_dist

Mapa de casos confirmados y activos por departamento

es_cant_casos_w_activo_dpto <-
  mapa_base_dpto + geom_sf(data = dpto_shape1,
                           aes(fill = cant_casos),
                           colour = "grey75") +
  scale_fill_fermenter(
    show.limits = TRUE,
    breaks = cant_casos_breaks,
    na.value = "grey60",
    palette = "YlOrRd",
    direction = 1,
    limits = c(1, max(dpto_shape1$cant_casos, na.rm = TRUE))
  ) +
  geom_point(
    data = datos_dptos,
    colour =  "Red",
    alpha = 0.5,
    aes(x = lon_dpto,
        y = lat_dpto,
        size = activo)
  ) +
  annotate(
    "rect",
    xmin = -63,
    xmax = -59.7,
    ymin = -26.1,
    ymax = -25.3,
    alpha = .2
  ) +
  annotate(
    "Text",
    x = -62.9,
    y = -25.7,
    label =  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)",
    size = 5,
    hjust = 0,
    family = "serif",
    color = "#4e4d47"
  ) +
  annotation_custom(
    bar_cant_casos_dpto,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  scale_size_continuous(
    name = "Total activos",
    range = c(2, 20),
    limits = c(1, max(datos_dptos$activo, na.rm = TRUE)),
    breaks = activo_breaks,
    labels = activo_labels
  ) +
   ## Para usar los brewer colors se implementa fill_fermenter o fill_distiller
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(fill = "Total confirmados",
       caption = fuentes,
       title = titullo_dpto,
       subtitle =  subtitulo)  +
  guides(
    fill = guide_bins(
      title.position = 'top',
      keywidth = unit(1.5, units = "cm"),
      keyheight = unit(0.8, units = "cm"),
      axis.linewidth = unit(0.4, units = "cm"),
      order = 1,
      show.limits = TRUE
    )
  ) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    panel.grid.major = element_line(size = 0.5, linetype = "solid", colour = "gray80"),
    plot.title = element_text(family = "serif",color = "#4e4d47", size = 20,  face = "bold"),
    plot.subtitle = element_text(family = "serif",color = "#4e4d47", size = 18),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = "right",
    legend.justification = c(0, 1),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "vertical",
  )
es_cant_casos_w_activo_dpto

Mapa de total de casos confirmados y activos por distrito

es_cant_casos_w_activo_dist <-
mapa_base_dpto + geom_sf(
  data = dist_shape1,
  aes(fill = cant_casos),
  colour = "grey75"
) +
  scale_fill_fermenter(
    show.limits = TRUE,
    breaks = cant_casos_breaks,
    na.value = "grey60",
    palette = "YlOrRd",
    direction = 1,
    limits = c(1, max(dist_shape1$cant_casos, na.rm = TRUE))
  )+
  geom_point(
    data = datos_dist,
    colour =  "Red",
    alpha = 0.5,
    aes(x = lon_dist,
        y = lat_dist,
        size = activo)
  ) +
  annotate("rect", xmin = -63, xmax =-59.7, ymin = -26.1, ymax = -25.3, alpha = .2) +
  annotate("Text", x=-62.9, y=-25.7, label=  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)", size=5,hjust = 0, family = "serif", color = "#4e4d47") +
 annotation_custom(
    bar_cant_casos_dist,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  )  +
  scale_size_continuous(
    name = "Total activos",
    range = c(2, 20),
    limits = c(1, max(datos_dist$activo, na.rm = TRUE)),
    breaks = activo_breaks,
    labels = activo_labels
  ) +
  ## Para usar los brewer colors se implementa fill_fermenter o fill_distiller
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(fill = "Total confirmados",
       caption = fuentes,
       title = titulo_dist,
       subtitle =  subtitulo)  +
  guides(
    fill = guide_bins(
      title.position = 'top',
      keywidth = unit(1.5, units = "cm"),
      keyheight = unit(0.8, units = "cm"),
      axis.linewidth = unit(0.4, units = "cm"),
      order = 1,
      show.limits = TRUE
    )
  ) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    panel.grid.major = element_line(size = 0.5, linetype = "solid", colour = "gray80"),
    plot.title = element_text(family = "serif",color = "#4e4d47", size = 20,  face = "bold"),
    plot.subtitle = element_text(family = "serif",color = "#4e4d47", size = 18),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = "right",
    legend.justification = c(0, 1),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "vertical",
  )
es_cant_casos_w_activo_dist
## Warning: Removed 34 rows containing missing values (geom_point).

Análisis de movilidad

Objetivos: Analizar la relación entre la variación de la movilidad y el varaciación de la cantidad de casos de covid en Py

Se analizaran las relaciones desde una ventana de días, entonces será necesario crear un dataframe en donde cada registro corresponda a una fecha y un departamento, con la cantidad de casos, y la variación de la movilidad a cada tipo de lugar

Mapas de movilidad de Paraguay por departamento

1- Mapa de variación de movilidad a comercios minoristas y lugares de recreación

inset_plot_tiendas_ocio <-# Plot linea de variacion de movilidad a comercios minoristas y lugares de recreación ----

  ggplotGrob(
    goog_py %>% filter(lugar == "Tiendas y espacio de ocio") %>%
      ggplot(aes(x = fecha, y = value)) + geom_line(alpha = 0.15, color = "gray65") + geom_smooth(
        se = FALSE,
        aes(
          x = fecha,
          y = value,
          color = lugar
        )
      )  +
      scale_y_continuous(
        name = "Variación porcentual de movilidad",
        breaks = c(-100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100),
        limits = c(-100, 100)
      ) + scale_x_date(
        name = "Fecha",
        date_breaks = "1 month",
        date_labels = "%b %d"
      ) + labs(
        color = "Tendencia de movilidad",
        title = "Variación del desplazamiento a Tiendas y espacio de ocio"
      ) +
      theme(
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "Gray75"),
        panel.grid.major.x = element_line(colour = "Gray75"),
        plot.title.position = "plot",
        plot.title = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 11,
          face = "bold",
          hjust = 0
        ),
        legend.position = c(0.3, 0.9),
        legend.text = element_text(
          color = "#4e4d47",
          family = "serif",
          hjust = 0,
          size = 8
        ),
        legend.title = element_blank(),
        legend.background = element_blank(),
        axis.title.x = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 7,
          face = "bold"
        ),
        axis.text.x = element_text(size = 7),
        axis.title.y = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 8,
          face = "bold"
        )) +
      annotate(
        geom = "text",
        x = fases_mov$fases_fechas,
        y = -99,
        label = fases_mov$fases_nombres,
        colour = "darkgray",
        angle = 90,
        hjust = 0,
        size = 2
      ) +
      geom_hline(yintercept = 0, colour = "red"))
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 351 rows containing non-finite values (stat_smooth).
es_var_tiendas_ocio <- # Mapa de variación de movilidad a comercios minoristas y lugares de recreación----
mapa_base_dpto + geom_sf(
  data = mov_py_actualizado,
  aes(fill = retail_and_recreation_percent_change_from_baseline),
  colour = "grey60"
) +
  scale_fill_fermenter(
    palette = "RdYlBu",
    breaks = c(-80, -60, -40, -20, 0, 20, 40, 60, 80),
    limits = c(-100, 100),
    direction = -1,
    oob = scales::squish
  ) +
   geom_label_repel(
    data = datos_dptos,
    aes(x = lon_dpto, y = lat_dpto, label = dpto_desc),
    fontface = "bold",
    size = 2,
    alpha = 0.5,
    segment.color = "black",
    force = 0,
    nudge_x = ifelse(
      datos_dptos$dpto_cod == "0",-0.4,
      ifelse(datos_dptos$dpto_cod == "11",-0.1, 0)
    ),
    nudge_y = 0
  ) +
  annotate(
    "rect",
    xmin = -63.1,
    xmax = -59.8,
    ymin = -27.2,
    ymax = -26.5,
    alpha = .2
  ) +
  annotate(
    "Text",
    x = -63,
    y = -26.8,
    label =  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)",
    size = 5,
    hjust = 0,
    family = "serif",
    color = "#4e4d47"
  ) +
annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  ) +
  annotation_custom(
    bar_cant_casos_dpto,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  annotation_custom(
    inset_plot_tiendas_ocio,
    xmin = -53.7,
    xmax = -50,
    ymin = -22,
    ymax = -18.8
  ) +
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(
    fill = "Cambio porcentual de desplazamiento",
    caption = fuentes,
    title = "Variación del desplazamiento a tiendas y espacios de ocio",
    subtitle =  subtitulo
  ) +
  guides(fill = guide_bins(
    title.position = 'top',
    keywidth = unit(1, units = "cm"),
    show.limits = TRUE
  )) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(
      size = 0.5,
      linetype = "solid",
      colour = "gray80"
    ),
    plot.background = element_blank(),
    plot.title = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 20,
      face = "bold"
    ),
    plot.subtitle = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 16
    ),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = c(0.02, 0.2),
    legend.justification = c(0, 0),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "horizontal",
  )
es_var_tiendas_ocio

2- Mapa de variación de movilidad a supermercados y farmacias

inset_plot_super_farmacia <-# Plot linea de variacion de movilidad a comercios minoristas y lugares de recreación ----

  ggplotGrob(
    goog_py %>% filter(lugar == "Supermercados y farmacias") %>%
      ggplot(aes(x = fecha, y = value)) + geom_line(alpha = 0.15, color = "gray65") + geom_smooth(
        se = FALSE,
        aes(
          x = fecha,
          y = value,
          color = lugar
        )
      )  +
      scale_y_continuous(
        name = "Variación porcentual de movilidad",
        breaks = c(-100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100),
        limits = c(-100, 100)
      ) + scale_x_date(
        name = "Fecha",
        date_breaks = "1 month",
        date_labels = "%b %d"
      ) + labs(
        color = "Tendencia de movilidad",
        title = "Variación del desplazamiento a supermercados y farmacias"
      ) +
      theme(
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "Gray75"),
        panel.grid.major.x = element_line(colour = "Gray75"),
        plot.title.position = "plot",
        plot.title = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 11,
          face = "bold",
          hjust = 0
        ),
        legend.position = c(0.3, 0.9),
        legend.text = element_text(
          color = "#4e4d47",
          family = "serif",
          hjust = 0,
          size = 8
        ),
        legend.background = element_blank(),
        legend.title = element_blank(),
        axis.title.x = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 7,
          face = "bold"
        ),
        axis.text.x = element_text(size = 7),
        axis.title.y = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 8,
          face = "bold"
        )) +
      annotate(
        geom = "text",
        x = fases_mov$fases_fechas,
        y = -99,
        label = fases_mov$fases_nombres,
        colour = "darkgray",
        angle = 90,
        hjust = 0,
        size = 2
      ) +
      geom_hline(yintercept = 0, colour = "red"))
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 509 rows containing non-finite values (stat_smooth).
es_var_super_farmacia <- # Mapa de variación de movilidad a comercios minoristas y lugares de recreación----
mapa_base_dpto + geom_sf(
  data = mov_py_actualizado,
  aes(fill = grocery_and_pharmacy_percent_change_from_baseline),
  colour = "grey60"
) +
  scale_fill_fermenter(
    palette = "RdYlBu",
    breaks = c(-80, -60, -40, -20, 0, 20, 40, 60, 80),
    limits = c(-100, 100),
    direction = -1,
    oob = scales::squish
  ) +
   geom_label_repel(
    data = datos_dptos,
    aes(x = lon_dpto, y = lat_dpto, label = dpto_desc),
    fontface = "bold",
    size = 2,
    alpha = 0.5,
    segment.color = "black",
    force = 0,
    nudge_x = ifelse(
      datos_dptos$dpto_cod == "0",-0.4,
      ifelse(datos_dptos$dpto_cod == "11",-0.1, 0)
    ),
    nudge_y = 0
  ) +
  annotate(
    "rect",
    xmin = -63.1,
    xmax = -59.8,
    ymin = -27.2,
    ymax = -26.5,
    alpha = .2
  ) +
  annotate(
    "Text",
    x = -63,
    y = -26.8,
    label =  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)",
    size = 5,
    hjust = 0,
    family = "serif",
    color = "#4e4d47"
  ) +
annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  ) +
  annotation_custom(
    bar_cant_casos_dpto,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  annotation_custom(
    inset_plot_super_farmacia,
    xmin = -53.7,
    xmax = -50,
    ymin = -22,
    ymax = -18.8
  ) +
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(
    fill = "Cambio porcentual de desplazamiento",
    caption = fuentes,
    title = "Variación del desplazamiento a supermercados y farmacias",
    subtitle =  subtitulo
  ) +
  guides(fill = guide_bins(
    title.position = 'top',
    keywidth = unit(1, units = "cm"),
    show.limits = TRUE
  )) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(
      size = 0.5,
      linetype = "solid",
      colour = "gray80"
    ),
    plot.background = element_blank(),
    plot.title = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 20,
      face = "bold"
    ),
    plot.subtitle = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 16
    ),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = c(0.02, 0.2),
    legend.justification = c(0, 0),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "horizontal",
  )
es_var_super_farmacia

3- Mapa de variación de movilidad a Parques

inset_plot_parque <-# Plot linea de variacion de movilidad a comercios minoristas y lugares de recreación ----

  ggplotGrob(
    goog_py %>% filter(lugar == "Parques") %>%
      ggplot(aes(x = fecha, y = value)) + geom_line(alpha = 0.15, color = "gray65") + geom_smooth(
        se = FALSE,
        aes(
          x = fecha,
          y = value,
          color = lugar
        )
      )  +
      scale_y_continuous(
        name = "Variación porcentual de movilidad",
        breaks = c(-100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100),
        limits = c(-100, 100)
      ) + scale_x_date(
        name = "Fecha",
        date_breaks = "1 month",
        date_labels = "%b %d"
      ) + labs(
        color = "Tendencia de movilidad",
        title = "Variación del desplazamiento a parques"
      ) +
      theme(
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "Gray75"),
        panel.grid.major.x = element_line(colour = "Gray75"),
        plot.title.position = "plot",
        plot.title = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 11,
          face = "bold",
          hjust = 0
        ),
        legend.background = element_blank(),
        legend.position = c(0.3, 0.9),
        legend.text = element_text(
          color = "#4e4d47",
          family = "serif",
          hjust = 0,
          size = 8
        ),
        legend.title = element_blank(),
        axis.title.x = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 7,
          face = "bold"
        ),
        axis.text.x = element_text(size = 7),
        axis.title.y = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 8,
          face = "bold"
        )) +
      annotate(
        geom = "text",
        x = fases_mov$fases_fechas,
        y = -99,
        label = fases_mov$fases_nombres,
        colour = "darkgray",
        angle = 90,
        hjust = 0,
        size = 2
      ) +
      geom_hline(yintercept = 0, colour = "red"))
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 887 rows containing non-finite values (stat_smooth).
es_var_mov_parque <- # Mapa de variación de movilidad a comercios minoristas y lugares de recreación----
mapa_base_dpto + geom_sf(
  data = mov_py_actualizado,
  aes(fill = parks_percent_change_from_baseline),
  colour = "grey60"
) +
  scale_fill_fermenter(
    palette = "RdYlBu",
    breaks = c(-80, -60, -40, -20, 0, 20, 40, 60, 80),
    limits = c(-100, 100),
    direction = -1,
    oob = scales::squish
  ) +
   geom_label_repel(
    data = datos_dptos,
    aes(x = lon_dpto, y = lat_dpto, label = dpto_desc),
    fontface = "bold",
    size = 2,
    alpha = 0.5,
    segment.color = "black",
    force = 0,
    nudge_x = ifelse(
      datos_dptos$dpto_cod == "0",-0.4,
      ifelse(datos_dptos$dpto_cod == "11",-0.1, 0)
    ),
    nudge_y = 0
  ) +
  annotate(
    "rect",
    xmin = -63.1,
    xmax = -59.8,
    ymin = -27.2,
    ymax = -26.5,
    alpha = .2
  ) +
  annotate(
    "Text",
    x = -63,
    y = -26.8,
    label =  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)",
    size = 5,
    hjust = 0,
    family = "serif",
    color = "#4e4d47"
  ) +
annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  ) +
  annotation_custom(
    bar_cant_casos_dpto,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  annotation_custom(
    inset_plot_parque,
    xmin = -53.7,
    xmax = -50,
    ymin = -22,
    ymax = -18.8
  ) +
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(
    fill = "Cambio porcentual de desplazamiento",
    caption = fuentes,
    title = "Variación del desplazamiento a parques",
    subtitle =  subtitulo
  ) +
  guides(fill = guide_bins(
    title.position = 'top',
    keywidth = unit(1, units = "cm"),
    show.limits = TRUE
  )) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(
      size = 0.5,
      linetype = "solid",
      colour = "gray80"
    ),
    plot.background = element_blank(),
    plot.title = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 20,
      face = "bold"
    ),
    plot.subtitle = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 16
    ),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = c(0.02, 0.2),
    legend.justification = c(0, 0),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "horizontal",
  )
es_var_mov_parque

4- Mapa de variación de movilidad a estaciones de transporte

inset_plot_mov_estaciones <-# Plot linea de variacion de movilidad a comercios minoristas y lugares de recreación ----

  ggplotGrob(
    goog_py %>% filter(lugar == "Estaciones de transporte") %>%
      ggplot(aes(x = fecha, y = value)) + geom_line(alpha = 0.15, color = "gray65") + geom_smooth(
        se = FALSE,
        aes(
          x = fecha,
          y = value,
          color = lugar
        )
      )  +
      scale_y_continuous(
        name = "Variación porcentual de movilidad",
        breaks = c(-100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100),
        limits = c(-100, 100)
      ) + scale_x_date(
        name = "Fecha",
        date_breaks = "1 month",
        date_labels = "%b %d"
      ) + labs(
        color = "Tendencia de movilidad",
        title = "Variación del desplazamiento a estaciones de transporte"
      ) +
      theme(
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "Gray75"),
        panel.grid.major.x = element_line(colour = "Gray75"),
        plot.title.position = "plot",
        plot.title = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 11,
          face = "bold",
          hjust = 0
        ),
        legend.background = element_blank(),
        legend.position = c(0.3, 0.9),
        legend.text = element_text(
          color = "#4e4d47",
          family = "serif",
          hjust = 0,
          size = 8
        ),
        legend.title = element_blank(),
        axis.title.x = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 7,
          face = "bold"
        ),
        axis.text.x = element_text(size = 7),
        axis.title.y = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 8,
          face = "bold"
        )) +
      annotate(
        geom = "text",
        x = fases_mov$fases_fechas,
        y = -99,
        label = fases_mov$fases_nombres,
        colour = "darkgray",
        angle = 90,
        hjust = 0,
        size = 2
      ) +
      geom_hline(yintercept = 0, colour = "red"))
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 1775 rows containing non-finite values (stat_smooth).
es_var_mov_estaciones <- # Mapa de variación de movilidad a comercios minoristas y lugares de recreación----
mapa_base_dpto + geom_sf(
  data = mov_py_actualizado,
  aes(fill = transit_stations_percent_change_from_baseline),
  colour = "grey60"
) +
  scale_fill_fermenter(
    palette = "RdYlBu",
    breaks = c(-80, -60, -40, -20, 0, 20, 40, 60, 80),
    limits = c(-100, 100),
    direction = -1,
    oob = scales::squish
  ) +
   geom_label_repel(
    data = datos_dptos,
    aes(x = lon_dpto, y = lat_dpto, label = dpto_desc),
    fontface = "bold",
    size = 2,
    alpha = 0.5,
    segment.color = "black",
    force = 0,
    nudge_x = ifelse(
      datos_dptos$dpto_cod == "0",-0.4,
      ifelse(datos_dptos$dpto_cod == "11",-0.1, 0)
    ),
    nudge_y = 0
  ) +
  annotate(
    "rect",
    xmin = -63.1,
    xmax = -59.8,
    ymin = -27.2,
    ymax = -26.5,
    alpha = .2
  ) +
  annotate(
    "Text",
    x = -63,
    y = -26.8,
    label =  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)",
    size = 5,
    hjust = 0,
    family = "serif",
    color = "#4e4d47"
  ) +
annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  ) +
  annotation_custom(
    bar_cant_casos_dpto,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  annotation_custom(
    inset_plot_mov_estaciones,
    xmin = -53.7,
    xmax = -50,
    ymin = -22,
    ymax = -18.8
  ) +
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(
    fill = "Cambio porcentual de desplazamiento",
    caption = fuentes,
    title = "Variación del desplazamiento a estaciones de transporte",
    subtitle =  subtitulo
  ) +
  guides(fill = guide_bins(
    title.position = 'top',
    keywidth = unit(1, units = "cm"),
    show.limits = TRUE
  )) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(
      size = 0.5,
      linetype = "solid",
      colour = "gray80"
    ),
    plot.background = element_blank(),
    plot.title = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 20,
      face = "bold"
    ),
    plot.subtitle = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 16
    ),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = c(0.02, 0.2),
    legend.justification = c(0, 0),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "horizontal",
  )
es_var_mov_estaciones

#### 5- Mapa de variación de movilidad a lugares de trabajo

inset_plot_mov_trabajo <-# Plot linea de variacion de movilidad a comercios minoristas y lugares de recreación ----

  ggplotGrob(
    goog_py %>% filter(lugar == "Lugares de trabajo") %>%
      ggplot(aes(x = fecha, y = value)) + geom_line(alpha = 0.15, color = "gray65") + geom_smooth(
        se = FALSE,
        aes(
          x = fecha,
          y = value,
          color = lugar
        )
      )  +
      scale_y_continuous(
        name = "Variación porcentual de movilidad",
        breaks = c(-100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100),
        limits = c(-100, 100)
      ) + scale_x_date(
        name = "Fecha",
        date_breaks = "1 month",
        date_labels = "%b %d"
      ) + labs(
        color = "Tendencia de movilidad",
        title = "Variación del desplazamiento a lugares de trabajo"
      ) +
      theme(
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "Gray75"),
        panel.grid.major.x = element_line(colour = "Gray75"),
        plot.title.position = "plot",
        plot.title = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 11,
          face = "bold",
          hjust = 0
        ),
        legend.background = element_blank(),
        legend.position = c(0.3, 0.9),
        legend.text = element_text(
          color = "#4e4d47",
          family = "serif",
          hjust = 0,
          size = 8
        ),
        legend.title = element_blank(),
        axis.title.x = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 7,
          face = "bold"
        ),
        axis.text.x = element_text(size = 7),
        axis.title.y = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 8,
          face = "bold"
        )) +
      annotate(
        geom = "text",
        x = fases_mov$fases_fechas,
        y = -99,
        label = fases_mov$fases_nombres,
        colour = "darkgray",
        angle = 90,
        hjust = 0,
        size = 2
      ) +
      geom_hline(yintercept = 0, colour = "red"))
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 107 rows containing non-finite values (stat_smooth).
es_var_mov_trabajo <- # Mapa de variación de movilidad a comercios minoristas y lugares de recreación----
mapa_base_dpto + geom_sf(
  data = mov_py_actualizado,
  aes(fill = workplaces_percent_change_from_baseline),
  colour = "grey60"
) +
  scale_fill_fermenter(
    palette = "RdYlBu",
    breaks = c(-80, -60, -40, -20, 0, 20, 40, 60, 80),
    limits = c(-100, 100),
    direction = -1,
    oob = scales::squish
  ) +
   geom_label_repel(
    data = datos_dptos,
    aes(x = lon_dpto, y = lat_dpto, label = dpto_desc),
    fontface = "bold",
    size = 2,
    alpha = 0.5,
    segment.color = "black",
    force = 0,
    nudge_x = ifelse(
      datos_dptos$dpto_cod == "0",-0.4,
      ifelse(datos_dptos$dpto_cod == "11",-0.1, 0)
    ),
    nudge_y = 0
  ) +
  annotate(
    "rect",
    xmin = -63.1,
    xmax = -59.8,
    ymin = -27.2,
    ymax = -26.5,
    alpha = .2
  ) +
  annotate(
    "Text",
    x = -63,
    y = -26.8,
    label =  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)",
    size = 5,
    hjust = 0,
    family = "serif",
    color = "#4e4d47"
  ) +
annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  ) +
  annotation_custom(
    bar_cant_casos_dpto,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  annotation_custom(
    inset_plot_mov_trabajo,
    xmin = -53.7,
    xmax = -50,
    ymin = -22,
    ymax = -18.8
  ) +
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(
    fill = "Cambio porcentual de desplazamiento",
    caption = fuentes,
    title = "Variación del desplazamiento a lugares de trabajo",
    subtitle =  subtitulo
  ) +
  guides(fill = guide_bins(
    title.position = 'top',
    keywidth = unit(1, units = "cm"),
    show.limits = TRUE
  )) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(
      size = 0.5,
      linetype = "solid",
      colour = "gray80"
    ),
    plot.background = element_blank(),
    plot.title = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 20,
      face = "bold"
    ),
    plot.subtitle = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 16
    ),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = c(0.02, 0.2),
    legend.justification = c(0, 0),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "horizontal",
  )
es_var_mov_trabajo

#### 6- Mapa de variación de permanencia en hogares

inset_plot_mov_hogares <-# Plot linea de variacion de movilidad a comercios minoristas y lugares de recreación ----

  ggplotGrob(
    goog_py %>% filter(lugar == "Permanencia en hogares") %>%
      ggplot(aes(x = fecha, y = value)) + geom_line(alpha = 0.15, color = "gray65") + geom_smooth(
        se = FALSE,
        aes(
          x = fecha,
          y = value,
          color = lugar
        )
      )  +
      scale_y_continuous(
        name = "Variación porcentual de movilidad",
        breaks = c(-100, -80, -60, -40, -20, 0, 20, 40, 60, 80, 100),
        limits = c(-100, 100)
      ) + scale_x_date(
        name = "Fecha",
        date_breaks = "1 month",
        date_labels = "%b %d"
      ) + labs(
        color = "Tendencia de movilidad",
        title = "Variación de la permanencia en hogares"
      ) +
      theme(
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "Gray75"),
        panel.grid.major.x = element_line(colour = "Gray75"),
        plot.title.position = "plot",
        plot.title = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 11,
          face = "bold",
          hjust = 0
        ),
        legend.background = element_blank(),
        legend.position = c(0.3, 0.9),
        legend.text = element_text(
          color = "#4e4d47",
          family = "serif",
          hjust = 0,
          size = 8
        ),
        legend.title = element_blank(),
        axis.title.x = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 7,
          face = "bold"
        ),
        axis.text.x = element_text(size = 7),
        axis.title.y = element_text(
          family = "serif",
          color = "#4e4d47",
          size = 8,
          face = "bold"
        )) +
      annotate(
        geom = "text",
        x = fases_mov$fases_fechas,
        y = -99,
        label = fases_mov$fases_nombres,
        colour = "darkgray",
        angle = 90,
        hjust = 0,
        size = 2
      ) +
      geom_hline(yintercept = 0, colour = "red"))
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 1748 rows containing non-finite values (stat_smooth).
## Warning: Removed 5 row(s) containing missing values (geom_path).
es_var_mov_hogares <- # Mapa de variación de movilidad a comercios minoristas y lugares de recreación----
mapa_base_dpto + geom_sf(
  data = mov_py_actualizado,
  aes(fill = residential_percent_change_from_baseline),
  colour = "grey60"
) +
  scale_fill_fermenter(
    palette = "RdYlBu",
    breaks = c(-80, -60, -40, -20, 0, 20, 40, 60, 80),
    limits = c(-100, 100),
    direction = -1,
    oob = scales::squish
  ) +
   geom_label_repel(
    data = datos_dptos,
    aes(x = lon_dpto, y = lat_dpto, label = dpto_desc),
    fontface = "bold",
    size = 2,
    alpha = 0.5,
    segment.color = "black",
    force = 0,
    nudge_x = ifelse(
      datos_dptos$dpto_cod == "0",-0.4,
      ifelse(datos_dptos$dpto_cod == "11",-0.1, 0)
    ),
    nudge_y = 0
  ) +
  annotate(
    "rect",
    xmin = -63.1,
    xmax = -59.8,
    ymin = -27.2,
    ymax = -26.5,
    alpha = .2
  ) +
  annotate(
    "Text",
    x = -63,
    y = -26.8,
    label =  "Obs: El análisis se realizó en base a \ncasos fuera de albergue (comunitarios)",
    size = 5,
    hjust = 0,
    family = "serif",
    color = "#4e4d47"
  ) +
annotation_custom(
    plot_curva_acu_acti,
    xmin = -53.7,
    xmax = -50,
    ymin = -24,
    ymax = -22
  ) +
  annotation_custom(
    bar_cant_casos_dpto,
    xmin = -53.7,
    xmax = -50,
    ymin = -28.2,
    ymax = -24
  ) +
  annotation_custom(
    inset_plot_mov_hogares,
    xmin = -53.7,
    xmax = -50,
    ymin = -22,
    ymax = -18.8
  ) +
  annotation_scale(location = "bl", width_hint = 0.4) +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    pad_x = unit(0.2, "in"),
    pad_y = unit(0.2, "in"),
    style = north_arrow_fancy_orienteering
  ) +
  labs(
    fill = "Cambio porcentual de desplazamiento",
    caption = fuentes,
    title = "Variación de permanencia en hogares",
    subtitle =  subtitulo
  ) +
  guides(fill = guide_bins(
    title.position = 'top',
    keywidth = unit(1, units = "cm"),
    show.limits = TRUE
  )) +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_line(
      size = 0.5,
      linetype = "solid",
      colour = "gray80"
    ),
    plot.background = element_blank(),
    plot.title = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 20,
      face = "bold"
    ),
    plot.subtitle = element_text(
      family = "serif",
      color = "#4e4d47",
      size = 16
    ),
    plot.caption = element_text(
      hjust = 0,
      size = 10,
      color = "#4e4d47",
      family = "serif"
    ),
    legend.title = element_text(
      face = "bold",
      color = "#4e4d47",
      family = "serif",
      hjust = 0,
      size = 16
    ),
    legend.position = c(0.02, 0.2),
    legend.justification = c(0, 0),
    legend.background = element_blank(),
    legend.text = element_text(
      color = "#4e4d47",
      family = "serif",
      size = 15
    ),
    legend.key = element_blank(),
    legend.direction = "horizontal",
  )
es_var_mov_hogares

### Análisis de datos de movilidad y COVID-19

Cual es el anális que se llevará acabo?

Primeramente hacemos algunas presunciones o preguntas

  • La variación de la cantidad de casos confirmados por día y/o por semana, en una determinada localidad se asocia directamente con la variación del movimiento de las personas en dicha localidad en ventanas temporales menores a 21 días?
  • Es posible medir el acatamiento de las restricciones gubernamentales en el marco de la COVID19 a través de la variación de la movilidad

Desarrollo:

Pregunta:¿Donde sería bueno medir esto?

Se tomaron Asunción, Central y Alto Paraná, para esta selección se tuvo en cuenta principalmente la siguiente limitación:

Los datos de movilidad para Paraguay están disponibles para departamentos, además no hay una continuidad en los datos de muchos de los departamentos debido principalmente a la densidad y derivado de eso la cantidad de datos También el tipo de actividad que se desarrolla en estos departamentos podría tomarse como similar teniendo en cuenta que cuentan con ciudades relativamente densas en el caso de Paraguay y que a raíz de eso existe predominancia del rubro comercial.

Entonces, para hacerlo comparable se debería evaluar:

1- Densidad poblacional 2- Rubro (Derivado de eso el impacto de las restricciones en la actividad de cada ciudad) 3- Se toman los datos a partir de la fecha en donde la cantidad de casos acumulados supera 50 registros comunitarios

El análisis se realizó de dos maneras: 1- Por departamentos: Se analizó la correlación existente entre la variación de la movilidad y la cantidad de casos confirmados por día

Análisis de correlación Asunción:

Matriz de correlación - Asunción

## Warning in write.csv2(asu_cov_casos, "asu_mob_exportar.csv", sep = ";"): attempt
## to set 'sep' ignored

###Análisis de correlación para el Departamento de Central

central_cov_casos <- datos_covid_py %>% filter(dpto_desc == "CENTRAL" & albergue == "NO")%>% group_by(fecha) %>% summarise("casesByDay" = n()) %>%  mutate(acum = cumsum(casesByDay), tres_antes =  rollsumr(casesByDay, k = 3, fill = NA), siete =rollsumr(casesByDay, k = 7, fill = NA), gr = log(tres_antes/3)/log(siete/7))
per_conf = 21

bar_casos_central <- # Grafico de barras por departamento para la cantidad de casos confirmados fuera de albergue----

    central_cov_casos %>%  ggplot(aes(
      x = fecha,
      y = casesByDay
    ))+ geom_bar(stat = "identity")
# Es necesario establecer una ventana de tiempo para el análisis y verificar que todas las variables de movilidad estén cuantificadas en cada fecha
# El periodo de análisis dependerá de los datos disponibles, sin en bargo sería bueno verificar que todas las fechas tengan las 6 categorías de movilidad 
central_mov <-
  goog_py %>% filter(location_code == "PY-11", #Filtro de asunción
                     date >= min(central_cov_casos$fecha - per_conf,# Filtro de fecha, todos los registros cuya fecha sea igual o posterior a 21 días antes de la fecha más antigua del set de datos de covid asuncion
                      date <= (max(central_cov_casos$fecha) - per_conf))) %>% mutate(fecha_ajustada = date + per_conf, fecha.Mov.Ori = date)# Se filtra para que el set de datos de movilidad no exceda 21 días antes de la última fecha de covid esto para hacerlos coincidir temporalmente

#Aqui ambos set de datos preparados sufriran una replicación pero bifurcada teniendo en cuenta cual será el análisis a llevar a cabo, de checho no tiene sentido evaluar los casos diarios con la movilidad de ese día?

#Primeramente recortaremos los datos de movilidad parar que coincidan con los datos de los casos, esto es que los datos entren dentro de la ventana de tiempo de casos de
#Luego llevaremos a cabo una union, en donde 
central_cor_set_shifted <- merge(central_mov,central_cov_casos,by.x = "fecha_ajustada", by.y = "fecha", all.x =  TRUE) %>% filter( acum != "NA") %>% pivot_wider(id_cols = c(gr,fecha.Mov.Ori,fecha_ajustada, casesByDay, acum),names_from = lugar, values_from = value) %>% mutate(mov_prom = (`Supermercados y farmacias` + `Lugares de trabajo` + Parques +`Estaciones de transporte`+ `Tiendas y espacio de ocio`)/5) %>% filter( gr != "NA" )


#CORRELACION  con fechas movidas
 cor_central <- central_cor_set_shifted %>% select(-fecha.Mov.Ori, -fecha_ajustada) %>% cor() %>% round(1)
 cor_central
##                             gr casesByDay acum Permanencia en hogares Parques
## gr                         1.0        0.0  0.0                   -0.1     0.1
## casesByDay                 0.0        1.0  0.9                   -0.2     0.3
## acum                       0.0        0.9  1.0                   -0.1     0.3
## Permanencia en hogares    -0.1       -0.2 -0.1                    1.0    -0.9
## Parques                    0.1        0.3  0.3                   -0.9     1.0
## Supermercados y farmacias  0.1        0.3  0.2                   -0.9     0.9
## Tiendas y espacio de ocio  0.1        0.3  0.2                   -0.9     1.0
## Lugares de trabajo         0.1        0.1  0.1                   -0.9     0.7
## Estaciones de transporte   0.1        0.1  0.1                   -0.9     0.9
## mov_prom                   0.1        0.2  0.2                   -1.0     0.9
##                           Supermercados y farmacias Tiendas y espacio de ocio
## gr                                              0.1                       0.1
## casesByDay                                      0.3                       0.3
## acum                                            0.2                       0.2
## Permanencia en hogares                         -0.9                      -0.9
## Parques                                         0.9                       1.0
## Supermercados y farmacias                       1.0                       1.0
## Tiendas y espacio de ocio                       1.0                       1.0
## Lugares de trabajo                              0.8                       0.8
## Estaciones de transporte                        0.9                       0.9
## mov_prom                                        1.0                       1.0
##                           Lugares de trabajo Estaciones de transporte mov_prom
## gr                                       0.1                      0.1      0.1
## casesByDay                               0.1                      0.1      0.2
## acum                                     0.1                      0.1      0.2
## Permanencia en hogares                  -0.9                     -0.9     -1.0
## Parques                                  0.7                      0.9      0.9
## Supermercados y farmacias                0.8                      0.9      1.0
## Tiendas y espacio de ocio                0.8                      0.9      1.0
## Lugares de trabajo                       1.0                      0.8      0.9
## Estaciones de transporte                 0.8                      1.0      0.9
## mov_prom                                 0.9                      0.9      1.0
central_correlation_jpg_byday <- ggcorrplot(cor_central, hc.order = TRUE, type = "lower",lab = TRUE)
central_correlation_jpg_byday

###Análisis de correlación para el Departamento de Alto Paraná

# Por semana----
 
altpar_cov_casos <- datos_covid_py %>% filter(dpto_desc == "ALTO PARANA" & albergue == "NO")%>% group_by(fecha) %>% summarise("casesByDay" = n()) %>%  mutate(acum = cumsum(casesByDay)) %>% filter(acum > 100)

# Visualizacion de casos por dia gr = (log(casesByDay/3)/log(casesByDay/7))
per_conf =21

# Es necesario establecer una ventana de tiempo para el análisis y verificar que todas las variables de movilidad estén cuantificadas en cada fecha
# El periodo de análisis dependerá de los datos disponibles, sin en bargo sería bueno verificar que todas las fechas tengan las 6 categorías de movilidad 
altpar_mov <-
  goog_py %>% filter(location_code == "PY-10", #Filtro de asunción
                     date >= min(altpar_cov_casos$fecha - per_conf,# Filtro de fecha, todos los registros cuya fecha sea igual o posterior a 21 días antes de la fecha más antigua del set de datos de covid asuncion
                      date <= (max(altpar_cov_casos$fecha) - per_conf))) %>% mutate(fecha_ajustada = date + per_conf, fecha.Mov.Ori = date)# Se filtra para que el set de datos de movilidad no exceda 21 días antes de la última fecha de covid esto para hacerlos coincidir temporalmente

#Aqui ambos set de datos preparados sufriran una replicación pero bifurcada teniendo en cuenta cual será el análisis a llevar a cabo, de checho no tiene sentido evaluar los casos diarios con la movilidad de ese día?

#Primeramente recortaremos los datos de movilidad parar que coincidan con los datos de los casos, esto es que los datos entren dentro de la ventana de tiempo de casos de
#Luego llevaremos a cabo una union, en donde 
altpar_cor_set_shifted <- merge(altpar_mov,altpar_cov_casos,by.x = "fecha_ajustada", by.y = "fecha", all.x =  TRUE) %>% filter( acum != "NA") %>% pivot_wider(id_cols = c(fecha.Mov.Ori,fecha_ajustada, casesByDay, acum),names_from = lugar, values_from = value) %>% mutate(mov_prom = (`Supermercados y farmacias` + `Lugares de trabajo` + Parques + `Tiendas y espacio de ocio`)/5) 


#CORRELACION  con fechas movidas
 cor_altpar <- altpar_cor_set_shifted %>% select(-fecha.Mov.Ori, -fecha_ajustada,-`Estaciones de transporte`) %>% cor() %>% round(2)
 cor_altpar
##                           casesByDay  acum Parques Tiendas y espacio de ocio
## casesByDay                      1.00 -0.02   -0.12                     -0.04
## acum                           -0.02  1.00    0.20                     -0.21
## Parques                        -0.12  0.20    1.00                      0.82
## Tiendas y espacio de ocio      -0.04 -0.21    0.82                      1.00
## Supermercados y farmacias       0.03 -0.18    0.61                      0.79
## Permanencia en hogares          0.13  0.22   -0.65                     -0.81
## Lugares de trabajo              0.04 -0.41   -0.14                      0.14
## mov_prom                       -0.04 -0.18    0.82                      0.95
##                           Supermercados y farmacias Permanencia en hogares
## casesByDay                                     0.03                   0.13
## acum                                          -0.18                   0.22
## Parques                                        0.61                  -0.65
## Tiendas y espacio de ocio                      0.79                  -0.81
## Supermercados y farmacias                      1.00                  -0.47
## Permanencia en hogares                        -0.47                   1.00
## Lugares de trabajo                            -0.05                  -0.50
## mov_prom                                       0.78                  -0.85
##                           Lugares de trabajo mov_prom
## casesByDay                              0.04    -0.04
## acum                                   -0.41    -0.18
## Parques                                -0.14     0.82
## Tiendas y espacio de ocio               0.14     0.95
## Supermercados y farmacias              -0.05     0.78
## Permanencia en hogares                 -0.50    -0.85
## Lugares de trabajo                      1.00     0.33
## mov_prom                                0.33     1.00
altpar_correlation_jpg_byday <- ggcorrplot(cor_altpar, hc.order = TRUE, type = "lower",lab = TRUE)
altpar_correlation_jpg_byday
## Exportando archivos

Exportando con ggsave, como se puede ver, el primer parámetro corresponde a el nombre que tendrá el archivo, es importante especificar la extensión, el segundo parámetro corres ponde al objeto, y luego para este caso se especifican el ancho y altura en pulgadas, vale la pena mencionar que es posible especificar la unidad de medida y también la calidad en dpi.

# # Exportando mapas interactivos en formato HTML----
#   saveWidget(mapa_cant_casos_fnegro_clabels_dpto, file="mapa_cant_casos_fnegro_clabels_dpto.html")  

ggsave( "es_cant_confirmados_dpto.png", es_cant_confirmados_dpto, width = 18, height = 11, dpi = 300)

ggsave( "es_cant_confirmados_dist.png", es_cant_confirmados_dist, width = 18, height = 11, dpi = 300)

ggsave( "es_cant_casos_w_activo_dpto.png", es_cant_casos_w_activo_dpto, width = 18, height = 11, dpi = 300)

ggsave( "es_cant_casos_w_activo_dist.png", es_cant_casos_w_activo_dist, width = 18, height = 11, dpi = 300)

#Movilidad
# Plots
  # Curvas de movilidad unificadas en el tiempo
  ggsave( "plot_py_mov_unificado.png", plot_py_mov_unificado, width = 18, height = 11, dpi = 300)
  
  ggsave( "plot_py_mov_individual.png", plot_py_mov_individual, width = 18, height = 11, dpi = 300)
  
# Mapas
  # Comercios minoristas
  ggsave( "es_var_tiendas_ocio.png", es_var_tiendas_ocio, width = 18, height = 11, dpi = 300)
  # Supermercados y farmacias
  ggsave( "es_var_super_farmacia.png", es_var_super_farmacia, width = 18, height = 11, dpi = 300)
  # Parques
  ggsave( "es_var_mov_parque.png", es_var_mov_parque, width = 18, height = 11, dpi = 300)
  # Estaciones de transporte
  ggsave( "es_var_mov_estaciones.png", es_var_mov_estaciones, width = 18, height = 11, dpi = 300)
  # Lugares de trabajo
  ggsave( "es_var_mov_trabajo.png", es_var_mov_trabajo, width = 18, height = 11, dpi = 300)
  # Permanencia en hogares
  ggsave( "es_var_mov_hogares.png", es_var_mov_hogares, width = 18, height = 11, dpi = 300)
# Graficos de correlacion
  #Asuncion
  ggsave( "asu_correlation_jpg_byday.png", asu_correlation_jpg_byday, dpi = 300)#asucion
  ggsave( "central_correlation_jpg_byday.png", central_correlation_jpg_byday, dpi = 300)#central
  ggsave( "altpar_correlation_jpg_byday.png", altpar_correlation_jpg_byday, dpi = 300)#central
  #Central

Bibliografía