library(introexercises)  # obtener datos para los ejericios 
library(learnr)          # create lessons from rmd 
library(gradethis)       # evaluate exercises
library(dplyr)           # wrangle data
library(flair)           # highlight code 
library(ggplot2)         # visualise data 
library(lubridate)       # work with dates
library(forcats)         # work with factors
library(fontawesome)     # for emojis 
library(scales)          # defining axes and units 
library(stringr)         # work with character strings
library(apyramid)        # creating demographic pyramids 
library(viridis)         # defining colour palettes 
library(janitor)         # clean data
library(tsibble)         # working with epiweeks
library(tidyr)           # clean data
# library(RMariaDB)        # connect to sql database

## set options for exercises and checking ---------------------------------------

## Define how exercises are evaluated 
gradethis::gradethis_setup(
  ## note: the below arguments are passed to learnr::tutorial_options
  ## set the maximum execution time limit in seconds
  exercise.timelimit = 60, 
  ## set how exercises should be checked (defaults to NULL - individually defined)
  # exercise.checker = gradethis::grade_learnr
  ## set whether to pre-evaluate exercises (so users see answers)
  exercise.eval = FALSE 
)

# ## event recorder ---------------------------------------------------------------
# ## see for details: 
# ## https://pkgs.rstudio.com/learnr/articles/publishing.html#events
# ## https://github.com/dtkaplan/submitr/blob/master/R/make_a_recorder.R
# 
# ## connect to your sql database
# sqldtbase <- dbConnect(RMariaDB::MariaDB(),
#                        user     = Sys.getenv("userid"),
#                        password = Sys.getenv("pwd"),
#                        dbname   = 'excersize_log',
#                        host     = "144.126.246.140")
# 
# 
# ## define a function to collect data 
# ## note that tutorial_id is defined in YAML
#     ## you could set the tutorial_version too (by specifying version:) but use package version instead 
# recorder_function <- function(tutorial_id, tutorial_version, user_id, event, data) {
#     
#   ## define a sql query 
#   ## first bracket defines variable names
#   ## values bracket defines what goes in each variable
#   event_log <- paste("INSERT INTO responses (
#                        tutorial_id, 
#                        tutorial_version, 
#                        date_time, 
#                        user_id, 
#                        event, 
#                        section,
#                        label, 
#                        question, 
#                        answer, 
#                        code, 
#                        correct)
#                        VALUES('", tutorial_id,  "', 
#                        '", tutorial_version, "', 
#                        '", format(Sys.time(), "%Y-%M%-%D %H:%M:%S %Z"), "',
#                        '", Sys.getenv("SHINYPROXY_PROXY_ID"), "',
#                        '", event, "',
#                        '", data$section, "',
#                        '", data$label,  "',
#                        '", paste0('"', data$question, '"'),  "',
#                        '", paste0('"', data$answer,   '"'),  "',
#                        '", paste0('"', data$code,     '"'),  "',
#                        '", data$correct, "')",
#                        sep = '')
# 
#     # Execute the query on the sqldtbase that we connected to above
#     rsInsert <- dbSendQuery(sqldtbase, event_log)
#   
# }
# 
# options(tutorial.event_recorder = recorder_function)
# hide non-exercise code chunks ------------------------------------------------
knitr::opts_chunk$set(echo = FALSE)
# data prep --------------------------------------------------------------------
combined <- rio::import(system.file("dat/linelist_combined_20141201.rds", package = "introexercises"))

combined <- combined %>% 
  mutate(week_onset      = yearweek(date_onset, week_start = 1), ## create week of onset variable  
         week_onset_date = as.Date(week_onset))                  ## create a date version 

# Create hospital weeks for heat plot
#####################################
# count cases by hospital-week
hospital_weeks_raw <- combined %>%
  count(hospital, week_onset)

# create longer dataset of all possible hospital-weeks
expanded <- hospital_weeks_raw %>%
  select(hospital, week_onset) %>%
  tidyr::expand(., week_onset, hospital)

# merge so that all hospital-weeks are represented in the data
hospital_weeks <- hospital_weeks_raw %>%
  right_join(expanded) %>%
  mutate(n = replace_na(n, 0)) %>%
  filter(week_onset >= ymd("2014-06-15"))




# 
# 
# count_data <- combined %>%
#   group_by(hospital, date_hospitalisation) %>%
#   summarize(n_cases = dplyr::n()) %>%
#   ungroup()
# # 
# # hospitalisation_week <- incidence(combined,
# #           date_index = date_hospitalisation,
# #           interval = "weeks")
# 
# demo_agg <- combined %>%
#   count(age_cat, , name = "cases") %>%
#   pivot_wider(
#     id_cols = age_cat,
#     names_from = ,
#     values_from = cases) %>%
#   rename(`missing_` = `NA`)
# 
# aggregate_data <- demo_agg %>%
#   pivot_longer(
#     col = c(female, male, missing_),            # cols to elongate
#     names_to = "",                # name for new col of categories
#     values_to = "counts") %>%           # name for new col of counts
#   mutate(
#      = na_if(, "missing_")) # convert "missing_" to NA

Introducción a R para Epidemiología Aplicada y Salud Pública

Bienvenido/a

Bienvenido al curso "Introducción a R para epidemiología aplicada", ofrecido por Epidemiología - una organización sin fines de lucro y el principal proveedor de formación, apoyo y herramientas de R para profesionales de primera línea de la salud pública.

knitr::include_graphics("images/logo.png", error = F)

Gráficos de calor,pirámides y curvas epidémicas

Este ejercicio se centra en gráficos más complejos de lo habitual en el trabajo de salud pública, como los gráficos de calor, las curvas epidémicas y las pirámides demográficas de edad/sexo.

Para ello vamos a utilizar conceptos que ya hemos revisado previamente, como:

Formato

Este ejercicio te guía a través de las tareas que debes realizar en RStudio en tu ordenador.

Obtener ayuda

Hay varias formas de obtener ayuda:

1) Busca la sección de ayuda (ver más abajo) 2) Pide ayuda al instructor/facilitador de tu curso en directo 3) Programa una llamada 1 a 1 con un instructor para "Tutoría del curso". 4) Publica una pregunta en La Comunidad de Applied Epi

Este es el aspecto que tendrá esa sección de ayuda:

r fontawesome::fa("lightbulb", fill = "gold") Haz clic para leer una pista

¡Aquí verás una pista útil!


r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

linelist %>% 
  filter(
    age > 25,
    district == "Bolo"
  )

Aquí tienes más explicaciones sobre por qué funciona la solución.


Responder a las preguntas del cuestionario te ayudará a comprender el material. Las respuestas no se registran.

Para practicar, responde a las siguientes preguntas:

quiz(
  question_radio("¿Cuándo debería mirar las secciones de ayuda en rojo?",
    answer("Tras intentar escribir el código yo mismo", correct = TRUE),
    answer("Antes de intentar escribir el código", correct = FALSE),
    correct = "Revisar el código de ejemplo después de intentar escribirlo usted mismo puede ayudarle a mejorar",
    incorrect = "Por favor, intente completar el ejercicio usted mismo, o use las pistas disponibles, antes de mirar la respuesta."
  )
)
question_numeric(
 "¿Qué tan ansioso estás por comenzar este tutorial? En una escala del 1 (menos ansioso) al 10 (más ansioso).?",
 answer(10, message = "Intenta no preocuparte, ¡te ayudaremos a conseguirlo!", correct = T),
 answer(9, message = "Intenta no preocuparte, ¡te ayudaremos a conseguirlo!", correct = T),
 answer(8, message = "Intenta no preocuparte, ¡te ayudaremos a conseguirlo!", correct = T),
 answer(7, message = "Intenta no preocuparte, ¡te ayudaremos a conseguirlo!", correct = T),
 answer(6, message = "De acuerdo, lo conseguiremos juntos", correct = T),
 answer(5, message = "De acuerdo, lo conseguiremos juntos", correct = T),
 answer(4, message = "¡Me gusta tu confianza!", correct = T),
 answer(3, message = "¡Me gusta tu confianza!", correct = T),
 answer(2, message = "¡Me gusta tu confianza!", correct = T),
 answer(1, message = "¡Me gusta tu confianza!", correct = T),
 allow_retry = TRUE,
 correct = "Gracias por compartir.",
 min = 1,
 max = 10,
 step = 1
)

Licencia

Por favor, envía un correo electrónico a contact@appliedepi.org si tienes preguntas sobre el uso de estos materiales.

Objetivos de aprendizaje

Durante este ejercicio se te pedirá que aprendas nuevas funciones leyendo su documentación. Esto es intencionado. Si te es confuso, pide ayuda a un instructor.

Refresca tu memoria

En primer lugar, refresquemos los conceptos básicos de ggplot2

quiz(
  question("¿Dónde se ubica la estética dinámica en el código ggplot?",
    allow_retry = TRUE,
    answer("dentro aes()", correct = T),
    answer("fuera aes()", message = "Las asignaciones dinámicas significan que la estética se asigna a una columna de los datos y la presentación puede variar para cada fila. Estas deben realizarse dentro de 'mapping = aes()'")
  ),
  question("¿La estética estática en la llamada inicial a ggplot() es heredada por geoms posteriores?",
    allow_retry = TRUE,
    answer("No", correct = TRUE),
    answer("Si", message = "No, por ejemplo, si establece el tamaño de los puntos en 3, las líneas tampoco tendrán el tamaño 3.")
  ),
  question("¿La estética dinámica en la llamada inicial a ggplot() es heredada por geoms posteriores?",
    allow_retry = TRUE,
    answer("No", message = "Sí, por ejemplo, la asignación de wt_kg al eje X es heredada por todas las geoms posteriores."),
    answer("Si", correct = TRUE)
  ),
    question("¿Deben realizarse ajustes al tema antes o después de configurar uno de los temas predeterminados?",
    allow_retry = TRUE,
    answer("Antes", message = "Un tema completo/predeterminado anulará todos los pequeños ajustes anteriores del tema."),
    answer("Después", correct = TRUE)
  ),
  question("¿Cuáles de los siguientes son temas prediseñados en ggplot?",
    allow_retry = TRUE,
    answer("theme_bw()", correct = TRUE),
    answer("theme_classic()", correct = TRUE),
    answer("theme_red()"),
    answer("scale_color_brewer()")
  ),
  question("¿Cómo esconderías la leyenda de un gráfico en ggplot?",
    allow_retry = TRUE,
    answer("theme(legend.title = 'element.blank()')"),
    answer("theme(legend.position = 'right')"),
    answer("theme(legend.position(none))"),
    answer("theme(legend.position = 'none')", correct = TRUE)
  ),
  question("¿Cómo configurarías tu leyenda para que aparezca en el centro de tu gráfico?",
    answer("theme(legend.position = 'middle')"),
    answer("theme(legend.position =  c(0.5,0.5))", correct = TRUE)
  )
)

Genial, ahora que hemos refrescado la memoria con el cuestionario anterior, empezaremos con nuestro primer tema, gráficos de calor.

Preparación

Abre tu proyecto ebola. Ejecuta todos los segmentos de código de tu script "ebola_sitrep.Rmd", para poder utilizarla base de datos combined.

Si no has completado el ejercicio de unión, o ves errores al intentar ejecutar tu script, puedes importar una "copia de seguridad". combined de la carpeta "data/clean/backup/" utilizando este comando:

combined <- import(here("data", "clean", "backup", "linelist_combined_20141201.rds"))

¡Habla con un instructor si tienes problemas para ejecutar tu script!

Paquetes

Añade los siguientes paquetes y cárgalos:

Recuerda colocar {tidyverse} como último paquete.

Tu sección de paquetes debería tener este aspecto

# Este comando para cargar los paquetes
pacman::p_load(
  rio,          # para importar datos
  here,         # para localizar archivos
  janitor,      # para limpieza de datos
  lubridate,    # para limpieza de dátiles
  epikit,       # para código en línea 
  escalas,      # para porcentajes y etiquetas de fecha
  officer,      # líneas fronterizas
  flextable,    # crear tablas HTML
  gtsummary,    # para mesas bonitas
  apiramid,     # para pirámides de edad/sexo
  RColorBrewer, # para paletas de colores
  gghighlight,  # resaltando partes de la trama
  ggExtra,      # funciones especiales de trazado
  viridis,      # para escalas de colores amigables para daltónicos
  tsibble,      # para epiweeks y otros análisis de series temporales
  tidyverse     # para gestión y visualización de datos
)

Nuevo segmento de código en tu R Markdown

Añade un nuevo segmento o sección "Gráficos de calor" cerca de la parte inferior de tu R Markdown.

# Gráfico de calor 

También puedes añadir un título en la parte de texto del R Markdown para delimitar este nuevo contenido.

Tu script "ebola_sitrep.Rmd" debería tener aproximadamente este esquema:

  1. YAML
  2. segmento de código de configuración
  3. Cargar paquetes
  4. Importar datos
  5. (Opcional) Análisis exploratorio
  6. Limpiar y exportar la lista de líneas de vigilancia
  7. Unir datos
  8. Resumen textual del brote
  9. Tablas resumen
  10. Gráficos
  11. Enfoque sobre un distrito concreto
  12. Transformar - tendencia de la evolución de los pacientes
  13. NUEVO - gráfico de calor
  14. (Opcional) area de zona de pruebas

Semanas Epidemiológicas o episemanas

Antes de continuar, merece la pena tratar la creación de "semanas" o "episemanas" (semanas epidemiológicas) en tu base de datos.

Para crear tablas de recuentos de casos por semanas, puede que necesites crear una nueva columna que contenga la "episemana", basada en una fecha como date_onset.

Hay muchas formas distintas de definir lo que es una "semana". Del mismo modo, algunos países formulan las episemanas a su criterio propio (por ejemplo, que empiecen los domingos, los lunes o los viernes). La mayor parte del mundo utiliza semanas que empiezan los lunes.

A continuación, te presentamos el enfoque más común y te proporcionamos una hoja de referencia para más detalles (en la carpeta "intro_course/learning materials").

Paquete {tsibble} para crear Semanas Epidémicas (por ejemplo, "2014 W44")

Recomendamos el paquete {tsibble} para crear valores de "Semana" epidémica. Este paquete se utiliza a menudo para el análisis de series temporales y tiene diversas funciones útiles sobre tiempo/fecha.

Su función yearweek() produce cómodamente semanas a partir de una columna de clase fecha. También puedes designar el inicio de semana con el argumento week_start =, donde 1 es para el lunes, 7 es para el domingo, etc.

Podemos crear una columna clase"yearweek" usando la función yearweek() dentro de mutate() llamada week_onset que mostrará valores como "2014 W44". Lo importante es que pueden utilizarse para dividir las filas en semanas discretas para una curva epidémica u otra tabla.

NO añadas este código a tu script, simplemente estamos mostrando cómo yearweek() altera el formato de las columnas.

# nueva columna para la semana de inicio de síntomas (como year-week)
combined %>% 
  mutate(week_onset = yearweek(date_onset, week_start = 1))

Los valores de week_onset tendrán este aspecto

combined %>% 
  mutate(week_onset = yearweek(date_onset, week_start = 1)) %>% 
  pull(week_onset) %>% 
  head()

Fechas de inicio de la semana epidémica

Para una máxima flexibilidad, ahora creamos una segunda columna que será una clase de fecha normal, que mostrará el fecha de inicio de la semana.

La creamos envolviendo el week_onset (clase "yearweek") en el la función del paquete {base} as.Date(). Las "semanas" se convertirán en fechas. Si seleccionaste week_start = 1 (lunes), la fecha devuelta será el lunes de esa semana (redondeando hacia abajo). Esto se guarda como una nueva columna, week_onset_date.

No añadas este código a tu script, es simplemente un ejemplo.

# new column for week of onset (as a date)
combined %>% 
  mutate(week_onset = yearweek(date_onset, week_start = 1)) %>% 
  mutate(week_onset_date = as.Date(week_onset)) 

Los valores de week_onset_date tendrán el siguiente aspecto:

combined %>% 
  mutate(week_onset = yearweek(date_onset, week_start = 1)) %>% 
  mutate(week_onset_date = as.Date(week_onset)) %>% 
  pull(week_onset_date) %>% 
  head()

De momento, añade este codigo anterior de mutate() al final de tu segmento de código "Unir", antes de combined sea exportado como archivo .rds. Más adelante, puedes optar por integrar estos comandos en la cadena de operadores pipe de limpieza de la base de datos de vigilancia.

Vuelve a ejecutar el comando "Unir" para que los cambios se apliquen en combined.

Estilo Avanzado de la función mutate()

mutate() permite más de una declaración dentro de sus paréntesis, separadas por comas. En un paso opcional, puedes insertar ambos cambios dentro del mismo mutate() comando.

combined <- combined %>% 
  mutate(week_onset      = yearweek(date_onset, week_start = 1), ## create week of onset variable  
         week_onset_date = as.Date(week_onset))                  ## create a date version 

¡Ahora tenemos dos versiones de columnas "epiweek" que podemos utilizar para hacer gráficos y tablas!

Gráficos de calor

Los gráficos de calor, también llamados "mosaicos de calor", son visualizaciones útiles cuando se trata de mostrar 3 variables (eje x, eje y y relleno).

Por ejemplo, tal vez quieras ver un desglose de cuántos casos se notificaron por semana en varios hospitales, para hacerte una idea de cómo ha progresado la epidemia en varios lugares. Veamos cómo funciona esto.

Supongamos que tenemos una base de datos llamada hospital_weeks que tiene las siguientes columnas (la vas a crear más adelante):

head(hospital_weeks)    # print the first 6 rows

Estos datos están en formato "largo", y cada posible semana-hospital tiene una fila.

Crea este base de datos en tu sección de código "Gráficos de calor" utilizando el código que aparece a continuación. No intentes entender este código (prueba nuestro curso avanzado de gestión de datos). En esencia, utiliza count() para crear una base de datos de recuentos de casos por semanas, luego crea una base de datos con todas las semanas-hospital posibles, y luego las une.

# conteo de casos por semanas-hospital
hospital_weeks_raw <- combined %>%
  count(hospital, week_onset_date) 

# crea una base de datos versión larga de todas posibles combinaciones semanas-hospital
expanded <- hospital_weeks_raw %>% 
  select(hospital, week_onset_date) %>%
  tidyr::expand(., week_onset_date, hospital)

# une todas las semanas-hospital disponibles en la base original con la base de todas las combinaciones de semanas-hospital
hospital_weeks <- hospital_weeks_raw %>%      
  right_join(expanded) %>%                            
  mutate(n = replace_na(n, 0)) %>% 
  filter(week_onset_date >= ymd("2014-06-15"))

Gráfico de calor básico

Revisar la documentación

r fontawesome::fa("eye", fill = "darkblue") Dedica unos minutos a revisar la documentación de la función geom_tile() que proviene de {ggplot2}. Puedes acceder a ella

1) Entrando en ?geom_tile en la Consola, o 2) Haciendo clic en el panel "Ayuda" (abajo a la derecha) e introduciendo geom_tile en la barra de búsqueda

Haz un gráfico de calor

Qué código ggplot() escribirías para crear un gráfico de calor que:

Añade un gráfico de calor a tu segmento de código "Gráficos de calor". Recuerda el formato de ggplot() en el que utilizamos + para añadir geom capas a la ggplot(). No olvides utilizar el argumento mapping = aes(x = , y = , fill = ) en tu ggplot() función.

r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡Pruébalo tú mismo primero!)

ggplot(data = hospital_weeks,
       mapping = aes(x = week_onset_date, y = hospital, fill = n)) +
  geom_tile()


quiz(
  question("¿Qué hospital informó la mayor cantidad de casos que comenzaron en septiembre?",
    allow_retry = TRUE,
    answer("SMMH"),
    answer("Hospital del Puerto", correct = T),
    answer("Hospital Militar"),
    answer("Hospital Central"),
    answer("Desconocido")
  ),
    question("¿Qué hospital informó la tasa de letalidad más alta?",
    allow_retry = TRUE,
    answer("SMMH", message = "No es posible conocer la tasa de letalidad de este gráfico"),
    answer("Hospital del Puerto", message = "No es posible conocer la tasa de letalidad de este gráfico"),
    answer("Hospital Militar", message = "No es posible conocer la tasa de letalidad de este gráfico"),
    answer("Hospital Central", message = "No es posible conocer la tasa de letalidad a partir de este gráfico"),
    answer("Desconocido", correct = T)
    )
)

Personalizar tu gráfico de calor

El gráfico es presentable, pero hagamos algunas cosas más para refinarlo.

Actualiza el código de tu gráfico de calor con lo siguiente:

1) Modifica hospital_weeks para crear un nuevo base de datos: hospital_weeks_adjusted. En este base de datos, asegúrate de que los valores NA de la columna hospital se recodifiquen como "No encontrado". y esta columna sea de clase "factor" con los niveles en este orden : "Hospital Central", "Hospital Militar", "Hospital del Puerto", "SMMH", "Otros", "Falta".

r fontawesome::fa("lightbulb", fill = "gold") Haz clic para leer una sugerencia

Dentro de un comando mutate(), utiliza el paquete {forcats} y su función fct_na_value_to_level() para redefinir la columna hospital para convertir NA en "No encontrado" (esto también convierte la columna en la clase "factor"). Dentro de la fct_na_value_to_level() incluye level = "No encontrado" para reasignar NA a "No encontrado".

Introduce esto en otro comando mutate, y utiliza la función fct_relevel() para ajustar el orden de los niveles.

Si necesitas un repaso sobre cómo cambiar NA los valores de las columnas de los factores y cómo reordenarlos, consulta el capítulo del Manual de R sobre factores.


r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)>

hospital_weeks_adjusted <- hospital_weeks %>% 
  mutate(hospital = fct_na_value_to_level(hospital, level = "No encontrado"),
         hospital = fct_relevel(
           hospital,
           "Hospital Central", "Hospital Militar", "Hospital del Puerto", "SMMH", "Otro", "No encontrado"))




2) En el ggplot de la nueva base de datos ajustada cambia el esquema de color a una escala de color divergente. Proporciona una escala que se extienda desde "azul cielo" en valores bajos hasta "tomate" en valores altos.

r fontawesome::fa("lightbulb", fill = "gold") Haz clic para leer una sugerencia

Añade scale_fill_gradient() a la ggplot() comando. Consulta la documentación de ayuda con ?scale_fill_gradient para entender los argumentos.


r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)>

ggplot(
  data = hospital_weeks_adjusted,
  mapping = aes(
    x = week_onset_date,
    y = hospital,
    fill = n)) +
  geom_tile()+
  scale_fill_gradient(low = "skyblue", high = "tomato")




3) Las marcas predeterminadas del eje x se muestran cada mes. Cámbielos para que se muestren de manera eficiente cada dos semanas.

r fontawesome::fa("lightbulb", fill = "gold") Haz clic para leer una sugerencia

Agregue scale_x_date() al final del ggplot y proporcione los argumentos date_breaks = "2 semanas" (por ejemplo) y labels = scales::label_date_short().


r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

ggplot(
  data = hospital_weeks_adjusted,
  mapping = aes(
    x = week_onset_date,
    y = hospital, fill = n)) +
  geom_tile()+
  scale_fill_gradient(low = "skyblue", high = "tomato")+
  scale_x_date(
    date_breaks = "2 weeks",
    labels = scales::label_date_short()
  )




4) Haz más presentables las etiquetas de los ejes y la leyenda.

r fontawesome::fa("lightbulb", fill = "gold") Haz clic para leer una sugerencia

Ajusta las etiquetas con labs(). Recuerda que el título de la leyenda en este caso se edita con el parámetro fill = porque fue creado por el relleno.


A continuación se muestra todo el código de la solución:

r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

# Edita el conjunto de datos
hospital_weeks_adjusted <- hospital_weeks %>% 
  mutate(hospital = fct_na_value_to_level(hospital, level = "No encontrado"),
         hospital = fct_relevel(
           hospital,
           "Hospital Central", "Hospital Militar", "Hospital del Puerto", "SMMH", "Otro", "No encontrado"))

# Gráfico con ajustes

ggplot(
  data = hospital_weeks_adjusted,
  mapping = aes(
    x = week_onset_date,
    y = hospital, fill = n)) +
  geom_tile()+
  scale_fill_gradient(low = "skyblue", high = "tomato")+
  scale_x_date(
    date_breaks = "2 weeks",
    labels = scales::label_date_short()
  )+
labs(
x = "Semana de inicio de los síntomas",
y = "Hospital",
fill = "Número de casos semanales")+
theme_minimal()


En el código de la solución, observe cómo podemos usar \n para dividir nuestras etiquetas en varias líneas. ¡Esta es una adición muy útil a nuestras etiquetas! Especialmente si queremos formatear los gráficos para que encajen bien en nuestro informe.

Configurar los datos, crear y personalizar un gráfico de calor

Tenga en cuenta que para hacer el gráfico de calor usando geom_tile(), tenía que proporcionar una base de datos con counts - datos agregados (no una lista de líneas). Hay formas de crear gráficos de calor a partir de datos de un listado (consulte las funciones geom_density()), pero no las cubriremos hoy.

Ahora que está familiarizado con la sintaxis de gráficos de mosaicos de calor, pasaremos a un ejercicio en el que creará una base de datos de conteos.

Intenta hacer lo siguiente por tu cuenta. Será un desafío, pero esta es una buena exposición a la programación R de nivel intermedio.

Primero, crea un nuevo fragmento de código para un resumen de edad y resultado, debajo de tu fragmento de mapas de calor. Haz lo siguiente:

r fontawesome::fa("check", fill = "red")Click to see a solution (try it yourself first!)

# crea una nueva base de datos haciendo un resumen de edad, sexo y fallecidos
age_outcome_summary <- combined %>%
  drop_na(sex, age_cat) %>%            # remueve valores faltantes
  group_by(age_cat, sex) %>%           # agrupa las filas por grupos de edad
  summarise(                           # comienza a crear un resumen por cada columna
    n = n(),
    n_death = sum(outcome == "Death", na.rm = TRUE), # total de filas con fallecidos
    pct_death = n_death / n)           # crea una columna con la proporción de fallecidos




2) Crea un gráfico de calor utilizando este nuevo base de datos. Haz que los mosaicos se rellenen según el porcentaje de resultados de Muerte, como se muestra en sex y age_cat.

Aquí tienes algunos retos adicionales, si tienes tiempo:

Añade esto a tu segmento de código "Resumen de resultados de edad",

r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

# hacer gráfico
ggplot(data = age_outcome_summary, # usar nuevo base de datos
       mapping  = aes(x = sex, # en el eje x
                     y = age_cat, # categoría de edad en el eje y
                     fill = pct_death)) + # el relleno (color) está sombreado por la proporción de muertos
  geom_tile() + # mostrar los datos como mosaicos  
  scale_fill_viridis( # ajustar la escala de colores
    option = "B", # elige cualquier opción
    limits = c(0, 1), # establece los límites para que vayan de 0 a 1
    labels = escalas::porcentaje, # ajusta los valores de la leyenda para que sean porcentajes
    direction = -1) + # invierte la dirección del gradiente de color

geom_text( # añade texto sobre los mosaicos
mapping  = aes(
label = scales::percent( # muestra porcentajes en lugar de proporciones
pct_death,
p = 0,1)))+
labs(x = "Sexo", # añadir etiquetas
y = "Categorías de edad",
fill = "Porcentaje de casos mortales")


quiz(
  question("¿Qué categoría combinada de edad y género tiene la proporción más alta de 'Muerte'?",
           answer("Hombres de 20-29", correct = TRUE),
           answer("Mujeres de 20-29"),
           answer("Hombres de 50-59"),
           answer("Mujeres de 40-49")
           ),
  question("Quien tiene las edades más altas?",
           answer("Mujeres"),
           answer("Hombres", correct = TRUE)
           )
)

Vale, ¡Buen trabajo! Has conseguido construir conjuntos de datos de resumen para un gráfico de calor y luego personalizar el gráfico para que puedas analizar la información.

Aunque no supieras cómo hacer todas esas ediciones al ggplot, el propósito de esa tarea era exponerte a un nuevo gráfico y a las posibilidades de personalización.

Siempre puedes encontrar código de ejemplo en el capítulo "Gráficos de calor" del manual Manual de Epi R y adquirir más práctica más adelante.

Curvas Epidemicas

Una curva epidémica (también conocida como "epicurva") es un gráfico epidemiológico básico que se suele utilizar para visualizar el patrón temporal de aparición de la enfermedad entre un conglomerado o epidemia de casos.

El análisis de la epicurva puede revelar tendencias temporales, valores atípicos, la magnitud del brote, el periodo de exposición más probable, intervalos de tiempo entre generaciones de casos, e incluso puede ayudar a identificar el modo de transmisión de una enfermedad no identificada (por ejemplo, fuente puntual, fuente común continua, propagación de persona a persona).

Vamos a construir curvas epidémicas con {ggplot2} que permite una personalización avanzada. Existe otra opción, la {incidence2} como se describe en el Manual Epi R que quizás sea más sencillo pero menos personalizable.

Como ya has hecho varios ejercicios utilizando {ggplot2} ya estás familiarizado con la sintaxis y la personalización que permite (temas, etiquetas, ejes, etc.).

Semanas o fechas

Cuando decidas hacer una curva epidémica, primero debes elegir entre si quieres mostrar la "semana" en el eje X (por ejemplo, "2014 W35") o algun formato de fecha del calendario.

Epicurvas con fechas del calendario

En su forma más básica, trazar curvas epidémicas en {ggplot2} consiste en utilizar la función ggplot geom_histogram() para mostrar la distribución de los valores de las fechas. Sin embargo, hay 3 cosas principales que debes tener en cuenta:

barras de tiempo simples

Crea una nueva sección y un segmento de código para "Epicurvas" debajo del "Resumen de resultados de edad".

Haz un histograma simple a partir del base de datos combined que muestre la distribución de date_onset utilizando intervalos de 5 días? El parámetro "binwidth" del histograma puede especificarse como un número (como una estética estática) para binwidth =. Además, limpia las etiquetas y utiliza el tema theme_bw().

Añade esto a tu segmento de código "Epicurvas". Recuerda que puedes utilizar en tu sección de código "Área de pruebas" para cualquier prueba de código que necesites hacer primero.

r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

ggplot(datos = combinados, mapeo = aes(x = fecha_inicio)) +
    geom_histogram(binwidth = 5) +
    labs(x = "Fecha de inicio", y = "Incidencia") +
    tema_bw()


Definir barras de tiempo manuales

Generar el gráfico utilizando la función bin_width() iniciará el primer intervalo en el primer caso. Como epidemiólogos y profesionales de la salud pública, normalmente necesitamos producir una curva epidémica que se ajuste a unas especificaciones muy concretas. Por ejemplo, una semana de 7 días que empiece en un día concreto. Del mismo modo, este enfoque no funciona para los intervalos mensuales, ¡Porque los meses pueden tener 28, 29, 30 o 31 días!

En lugar de utilizar binwidth = especificamos las barras del histograma con breaks =.

Un vector de fechas

breaks = espera un vector de fechas a utilizar como intervalos (puntos de corte).

Sabes que se puede hacer un vector utilizando la función c() por ejemplo c(54, 22, 89) o c("Paris", "Delhi", "Kigali").
Del mismo modo, podrías hacer un vector de fechas como c(ymd("2014-06-15"), ymd("2014-06-22"), ymd("2014-06-29")).
¡Pero sería un código muy largo para incluir todos los puntos de cortes semanales en este brote!

En su lugar, vamos a utilizar del paquete {base} la función seq.Date() para crear automáticamente un secuencia de fechas entre un from = valor, a un to = y en incrementos especificados a by =. Puedes especificar by valores como "día", "semana", "mes", etc.

La secuencia que se muestra a continuación utiliza las fechas del primer y último caso. Escribe este comando en tu fragmento de Prueba y ejecútalo para ver el resultado.

seq.Date(from = ymd("2014-05-06"),
         to = ymd("2014-11-28"),
         by = "week")

Esto se puede integrar en el comando de la epicurva. Pruébalo si quieres, pero que sepas que lo cambiaremos más adelante.

ggplot(data = combined, aes(x = date_onset))+
  geom_histogram(
    breaks = seq.Date(
      from = ymd("2014-05-06"),
      to = ymd("2014-11-28"),
      by = "week")) +
  labs(x = "Date of onset", y = "Incidence") +
  theme_bw()

Ahora estamos seguros de que las pausas utilizan intervalos de 7 días que empiezan y terminan en una fecha concreta. Pero normalmente necesitamos que los intervalos empiecen en Lunes (u otro día concreto).

Redondeo a la baja y a la alta en semanas

Para redondear fechas, te sugerimos que utilices floor_date() y su compañero ceiling_date() del paquete {lubridate}.

Prueba estos códigos en tu sección de Prueba. Simplemente muestran cómo funcionan estas funciones de forma aislada.

# Lunes antes del primer caso
floor_date(ymd("2014-05-06"), unit = "week", week_start = 1)

# Lunes después del último caso
ceiling_date(ymd("2014-11-28"), unit = "week", week_start = 1)

Al integrarlas en el seq.Date() podemos tener un vector con intervalos semanales a partir del lunes anterior al caso más antiguo.

Adapta tu seq.Date() para utilizar floor_date() y ceiling_date() como se indica a continuación.

# Secuenia de lunes antes del primer caso y después del último caso
seq.Date(from = floor_date(ymd("2014-05-06"), unit = "week", week_start = 1),
         to =   ceiling_date(ymd("2014-11-28"), unit = "week", week_start = 1),
         by =   "week")

Fechas de inicio y fin dinámicas

Por último, ¿Qué ocurre si nuestros datos se refrescan o actualizan con nuevos registros? ¡Esas fechas estáticas de nuestro código no serán exactas!

Sustituye las fechas estáticas por min() y max() utilizando la columna combined$date_onset. Porque este código no está dentro de una cadena de comandos (dplyr, operador pipe), tenemos que especificar la base de datos y la columna usando el nombre de la base de datos, el operador $ de índice y el nombre de la columna. No olvides que estas funciones necesitan na.rm = TRUE eliminar los valores que faltan.

seq.Date(from = floor_date(min(combined$date_onset, na.rm=T), unit = "week", week_start = 1),
         to =   ceiling_date(max(combined$date_onset, na.rm=T), unit = "week", week_start = 1),
         by =   "week")

De este modo, hemos conseguido que nuestro código sea más dinámico. El código dinámico garantiza que todos los gráficos y tablas se actualicen automáticamente.

Separa el comando pausas

La forma más fácil de organizar este comando es separar el comando seq.Date() del comando ggplot() comando.

Saca tu seq.Date() comando y muévelo antes del ggplot() guardándolo con un nombre como ebola_weeks que se proporciona al breaks = argumento en el geom_histogram().

# Define and save the vector
ebola_weeks <- seq.Date(
  from = floor_date(min(combined$date_onset, na.rm=T), unit = "week", week_start = 1),
  to =   ceiling_date(max(combined$date_onset, na.rm=T), unit = "week", week_start = 1),
  by =   "week")

Un último cambio antes de terminar es añadir closed = "left" a la función geom_histogram() que hace que las barras cuenten también los casos notificados en los días de interrupción.

# Define y guarda el vector con los intervalos de fecha
ebola_weeks <- seq.Date(
  from = floor_date(min(combined$date_onset, na.rm=T), unit = "week", week_start = 1),
  to =   ceiling_date(max(combined$date_onset, na.rm=T), unit = "week", week_start = 1),
  by =   "week")


# corre el codigo del grafico usando el vector

ggplot(data = combined, aes(x = date_onset)) +
  geom_histogram(
    breaks = ebola_weeks,
    closed = "left") +
  labs(x = "Date of onset", y = "Incidence") +
  theme_bw()

Compara este gráfico con el histograma realizado en el ejercicio anterior (7 días, sin pausas especificadas). ¿Ves alguna pequeña diferencia en cómo se han agrupado los casos? (mira al principio de la epidemia).

¡Lo anterior puede parecer mucho trabajo para producir sólo una lista de fechas! Pero date cuenta del control que tienes sobre el resultado. Y, ahora que tienes este código, será fácil modificarlo. Todo esto también está en el capítulo sobre Curvas Epidémicas del Manual de Epi R, ¡Para que te resulte fácil consultarlo!

Curva epidémica mensual

Ahora, crea una nueva curva epidémica en tu sección Curvas epidémicas, en la que los saltos de celda del histograma sean por mes desde el inicio del mes anterior al primer caso.

Una vez hecho esto ajusta el código para que el fill en las barras sea por hospital, y las etiquetas del eje x aparezcan cada mes.

r fontawesome::fa("lightbulb", fill = "gold") Haz clic para leer una sugerencia

Como queremos que los valores se agreguen por meses naturales, haz lo siguiente:

  • Establece el unit = a "mes" en los dos {lubridate} funciones y elimina el week_start argumentos (comprueba ?floor_date para más detalles).

  • Cambia el by = argumento de seq.Date() a "mes".

  • Cambia el nombre del vector ebola_weeks a ebola_months y cámbialo en ggplot() también.

  • Para fijar el relleno, pon fill = hospital en las asignaciones estéticas

  • Para ajustar las etiquetas del eje x, haz lo siguiente dentro de scale_x_date():

  • configurar labels = label_date_short()

  • configura date_breaks = "months"


r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)>

# define y guarda el vector (por meses)
ebola_meses <- seq.Fecha(
  from = floor_date(min(combined$date_onset, na.rm=T), unit = "mes"),
  to = ceiling_date(max(combined$date_onset, na.rm=T), unit = "mes"),
  por = "mes")

# ejecuta la trama con pausas mensuales, y rellena

ggplot(datos = combinados, mapeo = aes(x = fecha\_inicio, relleno = hospital)) +
geom\_histograma(pausas = ebola\_meses, cerrado = "izquierda") +
escala\_x\_fecha(
fecha\_rupturas = "meses",
etiquetas = escalas::etiqueta\_fecha\_corta())+
labs(x = "Fecha de inicio", y = "Incidencia") +
tema\_bw()


Mini gráficos (Facets)

Puede ser difícil leer esta epicurva y entender las tendencias de los casos de cada hospital.

Dividir los gráficos ("faceting" en inglés) significa crear "múltiples pequeños" (mini gráficos) para cada valor único en una columna (por ejemplo, hospital).

Crea estos mini-gráficos agregando facet_wrap() al ggplot con un +, y escribiendo una tilde ~ con el nombre de la columna que deseas utilizar. Por ejemplo:

# corre el código del gráfico con los intervalos mensuales y rellena por hospital las barras
ggplot(data = combined,
       mapping = aes(x = date_onset, fill = hospital)) +
  geom_histogram() +
  facet_wrap(~sex)

Ahora, crea facetas basadas en hospital para tu curva epidémica. Añade esto a tu R Markdown actualizando el código en tu fragmento de código "Epicurves.

r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

ggplot(datos = combinados, aes(x = inicio_fecha, relleno = hospital)) +
  geom_histogram(pausas = ebola_meses, cerrado = "izquierda") +
  escala_x_fecha(
    interrupciones_fecha = "meses",
    etiquetas = escalas::etiqueta_fecha_corta())+
  labs(x = "Fecha de inicio", y = "Incidencia") +
  tema_bw()+
  facet_wrap(~hospital)


Hay MUCHOS otros ajustes que puedes hacer a tu curva epidémica hecha con {ggplot2}. Muchos de ellos se describen en el Capítulo Curvas epidémicas del Manual de Epi.

Epicurvas con números de semana

Si quieres mostrar "números de semana", aprovecha el hecho de que creaste la columna week_onset al principio de este módulo. Aún así, habrá varios pasos:

1) Crea un vector de valores de la episemana (similar al que creamos en el apartado anterior, pero escrito como "2014 W35") 2) Crea un segundo vector de valores de fecha, a partir de los valores de la episemana 3) Utiliza count() en la columna discreta "año-semana" (por ejemplo week_onset) para obtener recuentos por semana 4) Genera un gráfico con geom_col() como se muestra a continuación

Añade otra sección de código dentro de la sección "Epicurvas" de tu script R Markdown para "Epicurva con números de semana".

A continuación, creamos un vector utilizando seq() - no seq.Date() porque estos valores son valores especiales de semana, no fechas.

## create vector of the unique weeks which occur
all_weeks <- seq(min(combined$week_onset, na.rm = TRUE),
                 max(combined$week_onset, na.rm = TRUE),
                 by = 1) 

all_weeks

Luego hacemos un segundo vector con estas semanas transformadas a fechas:

## change weeks to dates
all_weeks_date <-  as.Date(all_weeks)  

all_weeks_date

Luego utilizamos count() para obtener recuentos por episemana. La dirección as_tsibble() y fill_gaps(n = 0) son de la {tsibble} paquete y asegúrate de que las semanas sin casos estén representadas por 0, no descartadas.

counts <- combined %>%
  filter(!is.na(week_onset_date)) %>%
  count(week_onset_date) %>%
  as_tsibble() %>%     # ensure that all weeks are present in the table 
  fill_gaps(n = 0)     # filled in with 0 as necessary

Por último, trazamos los recuentos utilizando geom_col(). Observa:

## grafico 
counts %>%
  ggplot(aes(x = week_onset_date, y = n)) +        # n es el recuento 
  geom_col(width = 7) +                            # para columnas que dura la semana completa                             
  theme(axis.text.x = element_text(angle = 90)) +  # pone las etiquetas del eje x en un angulo de 90 grados
  scale_x_date(
    breaks = all_weeks_date,
    labels = all_weeks)

Pirámides de edad

Las pirámides de edades son una forma útil de ilustrar los datos demográficos, y pueden personalizarse de muchas formas distintas. Para nuestros fines, utilizaremos la función age_pyramid() de la base de datos apyramid paquete.

Crea un nuevo encabezado de sección "Pirámide de edades" y un segmento de código debajo de tu sección "Epicurvas".

{apyramid}

Es bueno practicar el aprendizaje de una función leyendo su documentación. La función age_pyramid() es bastante sencilla, ya que produce pirámides de edades de gran calidad con relativamente pocos argumentos.

Consulta la documentación de la función (?age_pyramid) y revisa los argumentos y ejemplos.

Ahora crea una pirámide de edades utilizando la función combined base de datos y el age_cat categorías, que se divide por sex. Añade esto a un nuevo fragmento de código "Pirámides de edad".

r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

pirámide_edad(datos = combinados,
            grupo_edad = gato_edad
            split_by = sexo)


También puedes representar los valores por proporción, en lugar de por recuento, e incluir una columna para los datos que faltan.

¿Puedes ahora, después de revisar la documentación, crear una pirámide de edad por sexo que incluya los datos que faltan (NA) y represente las proporciones?

r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

age_pyramid(datos = combined,
            age_group = age_cat,
            proportional  = TRUE
            na.rm = FALSE) 


Cuenta

En el ejemplo anterior, los datos están en formato de lista de líneas, donde cada fila es una observación única. Sin embargo, puedes recibir datos que ya estén agregados en recuentos.

Para practicar, crea este conjunto de datos demo_wide_counts ejecutando el siguiente comando en un nuevo segmento de código.

# Crea una base de datos de demostración
#############################
demo_wide_counts <- structure(list(
  age_cat = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, NA),
                      .Label = c("0-9", "10-19", "20-29", "30-39", "40-49", "50-59", "60-69", "70+"), class = "factor"),
  female = c(149L, 104L, 55L, 15L, 4L, 1L, NA, NA, NA),
  male = c(96L, 80L, 65L, 39L, 19L, 5L, 5L, 1L, NA),
  missing_sex = c(5L, 4L, 4L, 2L, NA, NA, 1L, NA, 26L)),
  row.names = c(NA, -9L),
  class = c("tbl_df", "tbl", "data.frame"))

Debería tener este aspecto

demo_wide_counts

Para ser compatible con age_pyramid() los datos deben estar en formato "largo", así que primero tendrías que transformar a formato largo para que la información esté en una sola columna (¡recuerda ordenar los datos!) y los recuentos estén en una sola columna.

Puedes utilizar este comando

demo_long_counts <- demo_wide_counts %>% 
  pivot_longer(
    col = c(female, male, missing_sex),    # cols a pasar a formato largo
    names_to = "sex",                      # nombre nuevo de las columnas de categorias
    values_to = "counts") %>%              # nombre nuevo de la columna de recuento
  mutate(
    sex = na_if(sex, "missing_sex"))       # convierte "missing_" to NA

Que da como resultado:

demo_long_counts

Ahora, utiliza lo que está escrito en la documentaciónage_pyramid() para ajustar el comando de trazado para que acepte recuentos.

r fontawesome::fa("check", fill = "red")Haz clic para ver una solución (¡pruébalo tú primero!)

pirámide_edad(datos = demo_cuentas_largas,
            grupo_edad = gato_edad, 
            split_by = sexo,
            recuento = recuentos)


Añadir ggplot a las pirámides de edad

Vamos a probarlo con nuestro combined objeto base de datos. Crea una pirámide de edad estratificada por sexo, añade esto al segmento de código "Pirámide de edad":

# Pirámides de edad con modificaciones de ggplot
age_pyramid(
     data = combined,
     age_group = age_cat,
     split_by = sex,
     proportional = TRUE,
     show_midpoint = FALSE)

También puedes añadir fácilmente comandos ggplot utilizando la función + a {apyramid} para añadir etiquetas y ajustar la escala de colores.

Pruébalo ahora:

# Pirámides de edad con modificaciones de ggplot
age_pyramid(
     data = combined,
     age_group = age_cat,
     split_by = sex,
     proportional = TRUE,
     show_midpoint = FALSE)+

  theme_minimal()+

  scale_fill_brewer(type = "qual", palette = 2)+

  labs(title = "Edad y sexo de los casos confirmados",
       y = "Proporción de todos los casos",
       x = "Grupo de edad",
       caption = "Brote de Ebola",
       fill = "Sexo")

Nota porque el ggplot() en la función {apyramid} age_pyramid() utiliza coord_flip() al añadir etiquetas en un gráfico age_pyramid(), tienes que proporcionar el x = y y = en el eje opuesto al que quieres que aparezcan.

Fin del ejercicio

¡Enhorabuena!

¡Has terminado el estudio del caso del ébola y has aprendido MUCHO código R durante estas sesiones!

Has aprendido mucho:

Asegúrate de guardar y montar tu script en un documento R Markdown.

Considera la posibilidad de añadir más contenido escrito entre los gráficos y las tablas, interpretando los resultados. Decide qué gráficos conservar o excluir. ¡La elección depende de ti!

Si quieres trabajar con un ejemplo de SIG con tu conjunto de datos sobre el ébola, pasa al Extra de más abajo...

Extra: Demostración SIG

Crear mapas en R

A continuación se muestra una breve demostración de cómo hacer mapas con R.

Puedes utilizar {ggplot2} para elaborar mapas y realizar análisis espaciales con relativa facilidad.
La lista de líneas Ébola incluye información sobre la latitud y la longitud de cada caso (lat y lon columnas). Hay archivos shape (capas de mapa) en la carpeta "ebola/data/shp". No vamos a dar una larga explicación de los SIG ni de los archivos shape, sin embargo, este ejercicio demostrará el trazado de los datos de la lista de líneas en mapas sencillos.

Si te interesa el SIG con R, pregunta a un instructor sobre nuestro curso avanzado.

Abre un nuevo script R. Puedes guardar el script R como "ebola_GIS_analyses.R" en tu proyecto "ebola" y en la subcarpeta "scripts".

#############################################
# GIS Ebola example
# Bonus section
# Your NAME here
#############################################

Primero tendrás que cargar paquetes específicos para el mapeo:

Añádelos a tu script de R, y recuerda incluir un encabezado de sección.

# cargar paquetes -----------------------------------------------------------
pacman::p_load(
  rio, 
  here,
  sf,           # para trabajar con datos geoespaciales
  ggspatial,    # para base de datos geoespacial y las flechas de norte
  raster,       # para formatear datos espaciales  
  tidyverse
)

A continuación, importa los datos necesarios. Utiliza el dataframe o base de datos combined. A continuación, importa los archivos shapefiles utilizando la función especial read_sf() especial.

# Importando los datos -------------------------------------------------------------

# importa la versión limpia y unida del listado de casos
combined <- import(here("data", "clean", "backup", "linelist_combined_20141201.rds"))

## importa el shapefile de los bordes del país 
shapefile <- read_sf(here("data", "shp", "sle_adm3.shp"))

Antes de hacer el gráfico tenemos que limpiar un poco los datos. Crea un objeto llamado districts para identificar los distintos valores de distrito presentes en el listado o base de datos. Esto nos permitirá filtrar los shapefiles para mantener sólo los distritos con casos presentes.

# Filtrar datos -------------------------------------------------------------

## distritos que nos interesan
districts <- combined %>%    
  distinct(admin3pcod) %>%   # identificar distintos valores presentes en la columna admin3pcod
  drop_na() %>%   # eliminar NA
  pull()          # saca los distintos valores a un vector

## filtrar shapefile para distritos de interés
shapefile <- shapefile %>% 
  filter(admin3Pcod %in% districts)

Ahora que has filtrado el shapefile a la región del brote, crea un ggplot básico utilizando los límites del admin.

# Gráfico sencillo de los polígonos del distrito. ---------------------------------------
## abrir un ggplot
shape_plot <- ggplot() + 

  ## agregue el archivo de forma en la parte superior
  geom_sf(data = shapefile, 
          fill = NA,         # sin relleno
          colour = "black")  # bordes negros
# presenta
shape_plot

Aquí asignamos al mapa el nombre shape_plot. Como lo asignamos a este nombre, no se presentará automáticamente en el panel del visor. Para ver el gráfico ejecuta el nombre del objeto. Esta es la diferencia entre guardar una salida (utilizando <- para crearla) y presentarla.

Este mapa se convertirá en una capa fundamental de nuestro gráfico. Al igual que podemos superponer geomas en ggplots que no son mapas, este shape_plot será la primera capa de nuestros gráficos de mapa.

Veremos 3 ejemplos sencillos de trazado:

Puntos

Para trazar puntos, elimina las filas con NAs presentes en lon y lat columnas. A continuación, define el sistema de referencia de coordenadas (SRC). No entraremos en detalles sobre esto aquí, pero ten en cuenta que varían según la región. Para más detalles, consulta el capítulo sobre SIG del Manual de Epi R.

############# PUNTOS  ##########################################################

combined_sf <- combined %>% 
  drop_na(lat, lon) %>% 
  st_as_sf(                                               
    # definir las coordenadas basadas en variables de latitud/longitud
    coords = c("lon", "lat"),                             
    # establecer el sistema de referencia de coordenadas en WGS84
    crs = 4326,                                           
    # no cambie las variables de cadena a factores 
    stringsAsFactors = FALSE                              
  )



# ver las primeras 10 filas, las primeras 5 columnas y la columna de geometría
combined_sf$geometry

Ten en cuenta que el objeto combined_sf creado ya no es una simple base de datos, sino un objeto espacial. En combined_sf hay un geometry que proporciona ubicaciones de casos con formato espacial. Ahora podemos trazar estos puntos en nuestra shape_plot utilizando la función geom_sf() (sf significa características espaciales).

## trazar puntos en las formas del distrito
shape_plot + 
  geom_sf(data = combined_sf)+
  labs(title = "Case locations")


## trazar puntos en las formas del distrito, coloreados según el resultado
shape_plot + 
     geom_sf(data = combined_sf,
             mapping = aes(color = fct_na_value_to_level(outcome)))+
     labs(color = "Condición de egreso",
          title = "Casos por localidad por condición de egreso") + 
  theme_minimal()

Muy interesante, ¿Verdad? Ya hemos creado nuestros primeros mapas de los datos de nuestra base de datos de casos de casos en R.

Observa que en el segundo gráfico, puedes colorear los casos por "Resultado", y también añadir un título. ¿Qué otras funciones podrías alterar utilizando la sintaxis ggplot que conoces?

Coropletas

A continuación añadimos una capa Coropleta para mostrar el recuento de casos de cada distrito afectado.

Para ello, debemos formatear nuestro objeto combined en una base de datos con los recuentos de casos por distrito. Utiliza columna admin3pcod para nuestros distritos. Una vez creado el objeto recuentos de casos, podemos fusionarlo con nuestro shapefile y luego trazar los recuentos utilizando la estética de relleno.

############# COROPLETAS ####################################################

## obtener recuentos de casos por distrito
case_counts <- combined %>%
  count(admin3pcod, name = "counts")



## agregar recuentos de casos al archivo de forma de los distritos 
shapefile <- left_join(shapefile, case_counts, by = c("admin3Pcod" = "admin3pcod"))

# mira el shapefile de los distritos y vea la nueva columna.
View(shapefile)

## grafica la capa de coropleta
ggplot() + 
  ## añade esta encima del shapefile
  geom_sf(data = shapefile, 
          # fill by case count
          aes(fill = counts),
          # black borders
          colour = "black")

Mapas base

Los mapas base proporcionan un mosaico subyacente con características geográficas como carreteras, ríos, bosques, etc. Debemos crear un "cuadro delimitador" basado en la extensión geográfica del shapefile. A continuación, se utiliza para recortar una fuente de acceso público como Google Maps u Open Street Map. Esto requiere una conexión a Internet.

############ MAPAS BASE #########################################################

# obtener el cuadro delimitador para el shapefile
bounding_box <- shapefile %>% 
  st_bbox()


# trazar un mapa base que incluya una barra de escala
basemap <- ggplot() +
# cambiar el cuadro delimitador a un objeto sf
# esto define el área para descargar mosaicos de mapas
  geom_sf(data = st_as_sfc(bounding_box)) +
  # descargar mosaicos de mapas y agregarlos al gráfico de mapa
  annotation_map_tile(
   # definir qué mosaicos de mapa usar
    type =  "cartolight",
  # definir carpeta para almacenar imágenes de mosaicos
    cachedir = here::here("data", "map_tiles"),
    # definir si se deben descargar mosaicos cada vez
    forcedownload = FALSE,
    # ocultar mensajes sobre el estado de descarga y el zoom
    progress = "none" )


# presenta el mapa base
basemap

Ahora podemos hacer un ggplot con los datos de la base de datos y graficarlos sobre el mapa base.

# grafica los casos encima del mapa base
basemap + 
  ## añade un shapefile encima
  geom_sf(data = shapefile, 
          # sin relleno
          fill = NA,
          # borders negros
          colour = "black") + 
  geom_sf(data = combined_sf,
          mapping = aes(color = fct_explicit_na(outcome)))+
  labs(color = "Condición de egreso",
       title = "Localidad de los casos, por condición de egreso")

¿Y si quisiéramos trazar las ubicaciones de los casos por resultado y mes? Utiliza la función facet_wrap() para formatear el mismo gráfico en mini-gráficos mensuales, de modo que podamos ver los casos a lo largo del tiempo.

# grafica los casos encima del mapa base
basemap + 
  ## añade el shapefile encima
  geom_sf(data = shapefile, 
          # sin relleno
          fill = NA,
          # bordes negros
          colour = "black") + 
  geom_sf(data = combined_sf,
          mapping = aes(color = fct_na_value_to_level(outcome)))+
  labs(color = "Condición de egreso",
       title = "Localidad de los casos, por condición de egreso")+
  facet_wrap(~month(date_onset, label = TRUE))

Los SIG con R son todo un tema y esto ha sido una simple demostración.

Dependiendo de la región, tendrás que descargarte shapefiles que representen tu región de interés.

¡No olvides guardar tu script!



appliedepi/introexercises documentation built on April 22, 2024, 1:01 a.m.