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
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)
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:
ggplot()
%>%
theme_minimal()
scale_y_continuous()
)Este ejercicio te guía a través de las tareas que debes realizar en RStudio en tu ordenador.
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 )
Por favor, envía un correo electrónico a contact@appliedepi.org si tienes preguntas sobre el uso de estos materiales.
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.
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.
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!
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 )
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:
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").
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()
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
.
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!
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"))
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
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) ) )
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".
NA
en "Falta" utilizando la función fct_na_value_to_level()
desde {forcats} dentro de mutate()
. Asegúrate de utilizar el parámetro level =
para cambiar NA
a No encontrado
entre comillas No encontrado
.fct_relevel()
también desde el paquete {forcats} en un mutate()
para proporcionar el orden de nivel deseador 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.
label_date_short()
de {scales}, que aprendió en el módulo escalas de ggplot. 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.
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:
combined
para crear una nueva base de datos llamada age_outcome_summary.outcome
) es igual a "Muerte" o "death".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:
"Percent of\ncases fatal"
scale_fill_viridis()
y elige un esquema de colores (por ejemplo, utiliza el argumento option = "B"
), y utiliza limits
argumento para establecer la escala limits = c(0, 1)
para alinearla con las proporciones.direction = -1
dentro de scale_fill_viridis()
para invertir la dirección de la paleta de colores, de modo que el porcentaje más alto sean colores oscuros y el más bajo sean colores claros.geom_text()
para mostrar el porcentaje en cada mosaico. Lee la documentación de geom_text
primero.labels = scales::percent
dentro de aes()
. Redondea a un decimal con el argumento accuracy = 0.1
.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.
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.).
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.
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:
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()
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 =
.
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).
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")
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.
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!
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()
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()
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.
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:
scale_x_date()
, los intervalos se definen como el vector de fechas porque es un eje de fechas## 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)
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".
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)
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)
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.
¡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...
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:
Para trazar puntos, elimina las filas con NA
s 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?
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")
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!
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.