Advertencias

Este documento está plegado de hipervínculos, ¡aprovéchalos!

En los recuadros con código en R, ¡usa el botón de la esquina superior derecha para copiar y pegar el contenido!

Los participantes deben conocer y comprometerse a cumplir el siguiete Código de Conducta. Tómate 5 minutos para leerlo

Objetivos

Analizar bases de datos provenientes de un sistema de vigilancia usando R
Generar gráficos con el paquete ggplot2 para describir en tiempo y espacio
Generar tablas con el paquete compareGroups para describir caracteristicas de persona
Estimar medidas de asociación con la función stats::glm para identificar factores de riesgo
Limpiar las salidas de la función stats::glm con el paquete epitidy
Construir una curva epidémica con el paquete incidence para calcular indicadores de tasa de crecimiento y tiempo de duplicación
Limpiar las salidas del paquete incidence con el paquete incidenceflow
Construir un canal endémico con epichannel para identificar zonas de éxito, seguridad, alerta y epidemia

INICIO

00 - ¿cómo empezar con R?

paso a paso

Desmitificar el uso de R creando gráficos y resolviendo ejercicios interactivos con Rstudio primers
Instalar R, Rstudio y un paquete con los siguientes pasos

Si llegas a presentar algún problema en tu computador, puedes acceder a través de tu buscador aquí

Instalar paquetes a usar aquí con la función install.packages(). (Pista: solo copiar y pegar!)

if(!require("remotes")) install.packages("remotes")
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("lubridate")) install.packages("lubridate")
if(!require("compareGroups")) install.packages("compareGroups")
if(!require("rio")) install.packages("rio")
if(!require("naniar")) install.packages("naniar")
if(!require("incidence")) install.packages("incidence")
if(!require("epiR")) install.packages("epiR")
if(!require("epiDisplay")) install.packages("epiDisplay")
if(!require("epitidy")) remotes::install_github("avallecam/epitidy")
if(!require("incidenceflow")) remotes::install_github("avallecam/incidenceflow")
if(!require("epichannel")) remotes::install_github("avallecam/epichannel")

Crear una carpeta de trabajo y proyecto en Rstudio. Más detalles del ¿qué, por qué y cómo? en video y diapositivas

miproyecto/
        |_ data-cruda/
        |_ data/
        |_ tabla/
        |_ archivos.R

algunos conceptos clave

Indentificar funciones, paquetes y datos en R. Ver una diapositiva resumen
Crear gráficos con la plantilla del paquete ggplot2. Ver una diapositiva resumen
Limpiar bases de datos usando verbos del paquete dplyr. Ver una diapositiva resumen
Usar el operador %>% llamado “pipe”. Ver un ejemplo aquí
Para más conceptos revisar el libro en línea The Epidemiologist R Handbook, específicamente el capítulo 3

01 - tiempo espacio persona

Material original del evento R aplicado a epidemiología

librerias a usar

library(tidyverse)
library(lubridate)
library(compareGroups)

importar base de datos

La fuente original de la base de datos es parte de un tutorial del RECON disponible aquí. (Pronto habrá un nuevo curso!)
Primero, importar la base de datos desde internet usando el paquete rio con la función import(), luego de una limpieza de datos

# https://github.com/avallecam/epiapli2019/blob/master/01-epi_descriptiva.R
ruta_limpio <- "https://github.com/avallecam/epiapli2019/raw/master/data/casoslimpio_20190916.rds"
casos_limpio <- rio::import(file = ruta_limpio)
casos_limpio

## # A tibble: 169 x 11
##    case_id generation date_of_infection date_of_onset date_of_hospitalisation
##    <fct>        <dbl> <date>            <date>        <date>                 
##  1 d1fafd           0 NA                2014-04-07    2014-04-17             
##  2 53371b           1 2014-04-09        2014-04-15    2014-04-20             
##  3 f5c3d8           1 2014-04-18        2014-04-21    2014-04-25             
##  4 6c286a           2 NA                2014-04-27    2014-04-27             
##  5 0f58c4           2 2014-04-22        2014-04-26    2014-04-29             
##  6 49731d           0 2014-03-19        2014-04-25    2014-05-02             
##  7 f9149b           3 NA                2014-05-03    2014-05-04             
##  8 881bd4           3 2014-04-26        2014-05-01    2014-05-05             
##  9 e66fa4           2 NA                2014-04-21    2014-05-06             
## 10 20b688           3 NA                2014-05-05    2014-05-06             
## # ... with 159 more rows, and 6 more variables: date_of_outcome <date>,
## #   outcome <fct>, gender <fct>, hospital <fct>, lon <dbl>, lat <dbl>

Segundo, identificar valores perdidos en base de datos usando el paquete naniar con la función miss_var_summary() y vis_miss()
Explorar las variables disponibles en bases de datos

casos_limpio %>% naniar::miss_var_summary()

## # A tibble: 11 x 3
##    variable                n_miss pct_miss
##    <chr>                    <int>    <dbl>
##  1 date_of_infection           63     37.3
##  2 date_of_outcome             56     33.1
##  3 case_id                      0      0  
##  4 generation                   0      0  
##  5 date_of_onset                0      0  
##  6 date_of_hospitalisation      0      0  
##  7 outcome                      0      0  
##  8 gender                       0      0  
##  9 hospital                     0      0  
## 10 lon                          0      0  
## 11 lat                          0      0

casos_limpio %>% naniar::vis_miss()

P1.1: ¿Por qué cree que hay valores perdidos en dichas variables?

distribución en tiempo

Tercero, crear un gráfico que permita ver la frecuencia de casos por unidad de tiempo.
Para ello, crear un histograma usando el paquete ggplot2 (tutorial primer)
Importante: Dentro de ggplot2 usamos + para conectar a sus componentes y capas, de una forma simmilar al %>%
yapa: ver video sobre ggplot2, diapositivas aquí

casos_limpio %>% 
  ggplot(aes(x = date_of_onset)) +
  geom_histogram(binwidth = 7, color="white") +
  scale_x_date(date_breaks = "7 day", date_labels = "%b-%d") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

P1.2: ¿Con qué otros análisis complementarías este gráfico?

distribución en espacio

Cuarto, evaluaremos si hay una región o lugar con una mayor cantidad de casos recientes.
Para ello, crear un diagrama de puntos. Pista: usar el cheatsheet para descubrir más geometrías!

casos_limpio %>% 
  ggplot(aes(x = lon, y = lat, colour = date_of_onset)) +
  geom_point() +
  scale_color_gradient(low = "red", high = "yellow", trans = "date") +
  theme_bw()

P1.3: ¿Crees que hay una agregación espacial de casos? ¿Qué más harías para rechazar dicha hipótesis?

tabla descriptiva

Quinto, describiremos las características de las persona en la base de datos
Seleccionar columnas outcome, gender, hospital con la función select() usando el paquete dplyr (tutorial primer)
Crear una tabla con formato usando el paquete compareGroups con la función compareGroups() y createTable()
Usar export2md() para ver en un formato en HTML
Transformar el orden de variables con la función mutate() y fct_infreq() usando el paquete dplyr (tutorial primer)

casos_limpio %>% 
  select(outcome, gender, hospital) %>% 
  mutate(hospital = fct_infreq(hospital),
         outcome = fct_infreq(outcome)) %>% 
  compareGroups(~., data = .) %>% 
  createTable() %>% 
  export2md()

Summary descriptives table
	[ALL]	N
	N=169
outcome:		169
NA	65 (38.5%)
Death	61 (36.1%)
Recover	43 (25.4%)
gender:		169
f	88 (52.1%)
m	81 (47.9%)
hospital:		169
NA	51 (30.2%)
Connaught Hospital	47 (27.8%)
other	25 (14.8%)
Military Hospital	24 (14.2%)
Princess Christian Maternity Hospital (PCMH)	12 (7.10%)
Rokupa Hospital	9 (5.33%)
Mitylira Hospital	1 (0.59%)

P1.4: ¿Preguntas?

02 - factores asociados

Material original del evento R aplicado a epidemiología y el paquete epitidy

librerias a usar

library(epitidy)

importar base de datos

Primero, inspeccionar documentación del estudio Whickham
Importar base de datos posterior a una limpieza de datos

# https://github.com/avallecam/epiapli2019/blob/master/02-clean_db.R
smoke_limpio <- "https://github.com/avallecam/epiapli2019/raw/master/data/smokeclean_20190906.rds"
smoke_clean <- rio::import(file = smoke_limpio)
smoke_clean

## # A tibble: 1,314 x 7
##    outcome smoker   age outcome_1 outcome_2 smoker_2 agegrp
##    <fct>   <fct>  <int>     <dbl> <fct>     <fct>    <fct> 
##  1 Alive   Yes       23         0 Alive     Yes      18-44 
##  2 Alive   Yes       18         0 Alive     Yes      18-44 
##  3 Dead    Yes       71         1 Dead      Yes      65+   
##  4 Alive   No        67         0 Alive     No       65+   
##  5 Alive   No        64         0 Alive     No       45-64 
##  6 Alive   Yes       38         0 Alive     Yes      18-44 
##  7 Alive   Yes       45         0 Alive     Yes      45-64 
##  8 Dead    No        76         1 Dead      No       65+   
##  9 Alive   No        28         0 Alive     No       18-44 
## 10 Alive   No        27         0 Alive     No       18-44 
## # ... with 1,304 more rows

tabla descriptiva

Segundo, crear una tabla con formato usando el paquete compareGroups
según la variable dependiente outcome, describir la distribución de tres variables smoker+agegrp+age
en la función compareGroups() especificar el argumento byrow = T
en la función createTable() especificar los argumentos show.all = T, sd.type = 2

tabla1 <- smoke_clean %>% 
  compareGroups(outcome~smoker+agegrp+age, data = ., byrow = T) %>% 
  createTable(show.all = T, sd.type = 2)

tabla1

## 
## --------Summary descriptives table by 'outcome'---------
## 
## _______________________________________________________ 
##              [ALL]       Alive       Dead     p.overall 
##             N=1314       N=945       N=369              
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## smoker:                                         0.003   
##     No    732 (55.7%) 502 (68.6%) 230 (31.4%)           
##     Yes   582 (44.3%) 443 (76.1%) 139 (23.9%)           
## agegrp:                                        <0.001   
##     18-44 624 (47.5%) 597 (95.7%) 27 (4.33%)            
##     45-64 447 (34.0%) 314 (70.2%) 133 (29.8%)           
##     65+   243 (18.5%) 34 (14.0%)  209 (86.0%)           
## age        46.9±17.4   40.1±13.8   64.4±12.8   <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

P2.1: ¿Qué tendencia logras identificar en las variables?

Tercero, exportar en formato XLSX con la función export2xls()

tabla1 %>% 
  export2xls("tabla/tabla02.xls")

medidas de asociación

Cuarto, estimar el riesgo de fallecimiento según el hábito de fumar.
A un diseño de estudio epidemiológico de cohorte le corresponde estimar el Riesgo Relativo (RR)
Para ello, construir modelos de regresión lineal con la función glm() con familia de distribución Poisson función de enlace log
Luego, limpiar salidas de modelos usando el paquete epitidy con la función epi_tidymodel_rr()

modelo simple

# smoke_clean %>% glimpse()

#simple
simple_model <- glm(outcome_1 ~ smoker, 
                    data = smoke_clean, 
                    family = poisson(link = "log"))
epi_tidymodel_rr(simple_model)

## # A tibble: 2 x 7
##   term        log.rr     se    rr conf.low conf.high p.value
##   <chr>        <dbl>  <dbl> <dbl>    <dbl>     <dbl>   <dbl>
## 1 (Intercept) -1.16  0.0659 0.314    0.275     0.357  0     
## 2 smokerYes   -0.274 0.107  0.760    0.615     0.937  0.0107

modelo múltiple

#multiple: controlar por confusión
multiple_model <- glm(outcome_1 ~ smoker + age, 
                      data = smoke_clean, 
                      family = poisson(link = "log"))
epi_tidymodel_rr(multiple_model)

## # A tibble: 3 x 7
##   term         log.rr      se      rr conf.low conf.high p.value
##   <chr>         <dbl>   <dbl>   <dbl>    <dbl>     <dbl>   <dbl>
## 1 (Intercept) -4.96   0.261   0.00702  0.00416    0.0116   0    
## 2 smokerYes    0.166  0.112   1.18     0.947      1.47     0.136
## 3 age          0.0644 0.00374 1.07     1.06       1.07     0

P2.2: ¿Cómo interpretas el cambio de RR entre muerte~fumar del modelo simple al modelo múltiple ajustado por edad?

exportar en formato XLSX usando el paquete writexl con la función write_xlsx()

epi_tidymodel_rr(multiple_model) %>% 
  writexl::write_xlsx("tabla/tabla03.xlsx")

métodos alternativo

# medidas de asociación ---------------------------------------------------

#epiR
library(epiR)
smoke_tabla1 <- with(smoke_clean,table(smoker_2,outcome_2)) %>% print()
epi.2by2(smoke_tabla1,method = "cohort.count")

#epiDisplay
library(epiDisplay)
smoke_tab2 <- with(smoke_clean,table(outcome,smoker)) %>% print()
cs(cctable = smoke_tab2)

# controlar por confusión -------------------------------------------------

#Mantel-Haenszel
smoke_tab3 <- with(smoke_clean,table(smoker_2,outcome_2,agegrp)) %>% print()
epi.2by2(smoke_tab3,method = "cohort.count")
mhor(mhtable=smoke_tab3,graph = F,design = "cohort")

P2.3: ¿Preguntas?

INTERMEDIO

RETORNAMOS EN 5 MINUTOS

03 - curva epidémica

Material ilustrativo disponible en la presentación Análisis de múltiples epidemias y prevalencias con R y purrr y extraido del paquete incidenceflow

librerias a usar

library(outbreaks) #sample data
library(incidence) #core functions
library(incidenceflow)

importar base de datos

Primero, acceder a la base de datos linelist dentro de la lista ebola_sim

ebola_sim$linelist %>% as_tibble()

## # A tibble: 5,888 x 11
##    case_id generation date_of_infection date_of_onset date_of_hospitalisation
##    <chr>        <int> <date>            <date>        <date>                 
##  1 d1fafd           0 NA                2014-04-07    2014-04-17             
##  2 53371b           1 2014-04-09        2014-04-15    2014-04-20             
##  3 f5c3d8           1 2014-04-18        2014-04-21    2014-04-25             
##  4 6c286a           2 NA                2014-04-27    2014-04-27             
##  5 0f58c4           2 2014-04-22        2014-04-26    2014-04-29             
##  6 49731d           0 2014-03-19        2014-04-25    2014-05-02             
##  7 f9149b           3 NA                2014-05-03    2014-05-04             
##  8 881bd4           3 2014-04-26        2014-05-01    2014-05-05             
##  9 e66fa4           2 NA                2014-04-21    2014-05-06             
## 10 20b688           3 NA                2014-05-05    2014-05-06             
## # ... with 5,878 more rows, and 6 more variables: date_of_outcome <date>,
## #   outcome <fct>, gender <fct>, hospital <fct>, lon <dbl>, lat <dbl>

A partir de la base de datos linelist, extraer el vector date_of_onset

dat <- ebola_sim$linelist$date_of_onset
enframe(dat,value = "date_of_onset") %>% select(date_of_onset)

## # A tibble: 5,888 x 1
##    date_of_onset
##    <date>       
##  1 2014-04-07   
##  2 2014-04-15   
##  3 2014-04-21   
##  4 2014-04-27   
##  5 2014-04-26   
##  6 2014-04-25   
##  7 2014-05-03   
##  8 2014-05-01   
##  9 2014-04-21   
## 10 2014-05-05   
## # ... with 5,878 more rows

Luego, configurar el vector aislado con la función incidence() y el argumento interval=7

i.7 <- incidence(dat, interval=7)
i.7

## <incidence object>
## [5888 cases from days 2014-04-07 to 2015-04-27]
## [5888 cases from ISO weeks 2014-W15 to 2015-W18]
## 
## $counts: matrix with 56 rows and 1 columns
## $n: 5888 cases in total
## $dates: 56 dates marking the left-side of bins
## $interval: 7 days
## $timespan: 386 days
## $cumulative: FALSE

Generar un gráfico con la función plot()

plot(i.7)

curva ascendente

Segundo, restringir los primeros 20 días con la expresión .[1:20]
Calcular tasa de crecimiento y tiempo de duplicación con la función fit()
Revisar las diapositivas 6 a 9 para interepretar los indicadores

f1 <- fit(i.7[1:20])
f1

## <incidence_fit object>
## 
## $model: regression of log-incidence over time
## 
## $info: list containing the following items:
##   $r (daily growth rate):
## [1] 0.03175771
## 
##   $r.conf (confidence interval):
##           2.5 %     97.5 %
## [1,] 0.02596229 0.03755314
## 
##   $doubling (doubling time in days):
## [1] 21.8261
## 
##   $doubling.conf (confidence interval):
##         2.5 %   97.5 %
## [1,] 18.45777 26.69823
## 
##   $pred: data.frame of incidence predictions (20 rows, 5 columns)

curva descendente

Tercero, identificar punto de corte de la curva con la función fit_optim_split()

f2 <- fit_optim_split(i.7)
f2

## $df
##         dates   mean.R2
## 1  2014-08-04 0.7650406
## 2  2014-08-11 0.8203351
## 3  2014-08-18 0.8598316
## 4  2014-08-25 0.8882682
## 5  2014-09-01 0.9120857
## 6  2014-09-08 0.9246023
## 7  2014-09-15 0.9338797
## 8  2014-09-22 0.9339813
## 9  2014-09-29 0.9333246
## 10 2014-10-06 0.9291131
## 11 2014-10-13 0.9232523
## 12 2014-10-20 0.9160439
## 13 2014-10-27 0.9071665
## 
## $split
## [1] "2014-09-22"
## 
## $fit
## <list of incidence_fit objects>
## 
## attr(x, 'locations'): list of vectors with the locations of each incidence_fit object
## 
## 'before'
## 'after'
## 
## $model: regression of log-incidence over time
## 
## $info: list containing the following items:
##   $r (daily growth rate):
##      before       after 
##  0.02982209 -0.01016191 
## 
##   $r.conf (confidence interval):
##              2.5 %       97.5 %
## before  0.02608945  0.033554736
## after  -0.01102526 -0.009298561
## 
##   $doubling (doubling time in days):
##   before 
## 23.24274 
## 
##   $doubling.conf (confidence interval):
##           2.5 %  97.5 %
## before 20.65721 26.5681
## 
##   $halving (halving time in days):
##    after 
## 68.21031 
## 
##   $halving.conf (confidence interval):
##          2.5 %   97.5 %
## after 62.86899 74.54349
## 
##   $pred: data.frame of incidence predictions (57 rows, 6 columns)
## 
## $plot

generar tablas resumen

Cuarta, transformar salidas a formato base de datos usando el paquete incidenceflow

curva ascendente

Con la función tidy_incidence() generamos una tabla solo con la tasa de crecimiento y tiempo de duplicación

f1 %>% tidy_incidence()

## # A tibble: 2 x 5
##   mark  parameter estimate conf_lower conf_upper
##   <chr> <chr>        <dbl>      <dbl>      <dbl>
## 1 1     r           0.0318     0.0260     0.0376
## 2 1     doubling   21.8       18.5       26.7

Con la función glance_incidence() generamos una tabla solo con las medidas de bondad de ajuste del modelo lineal

f1 %>% glance_incidence()

## # A tibble: 1 x 12
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.880         0.874 0.498      133. 9.81e-10     1  -13.4  32.8  35.7
## # ... with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

curva descendente

f2 %>% pluck("fit") %>% tidy_incidence()

## # A tibble: 4 x 5
##   mark   parameter estimate conf_lower conf_upper
##   <chr>  <chr>        <dbl>      <dbl>      <dbl>
## 1 before r           0.0298     0.0261    0.0336 
## 2 after  r          -0.0102    -0.0110   -0.00930
## 3 before doubling   23.2       20.7      26.6    
## 4 after  halving    68.2       62.9      74.5

f2 %>% pluck("fit") %>% glance_incidence()

## # A tibble: 2 x 13
##   names  r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC
##   <chr>      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl>
## 1 before     0.922         0.919 0.455      273. 2.95e-14     1  -14.8  35.5
## 2 after      0.951         0.949 0.155      578. 3.72e-21     1   15.4 -24.8
## # ... with 4 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>,
## #   nobs <int>

P3.1: ¿La tasa de crecimiento estimada antes y despúes del punto de corte es concordante con lo observado en el gráfico?

P3.2: ¿Preguntas?

04 - canal endémico

Material extraído del paquete epichannel

librerias a usar

library(epichannel)

importar base de datos

Primero, leer datos de vigilancia usando el paquete readr con la función read_csv()
Transformar algunas variables al formato adecuado

datos de enfermedad

denv <-
  readr::read_csv("https://predict.cdc.gov/api/v1/attachments/dengue%20forecasting%20project/iquitos_training_data.csv") %>%
  mutate(year = lubridate::year(week_start_date),
         epiweek = lubridate::epiweek(week_start_date)) %>%
  mutate(adm="iquitos") %>%
  # cases per season - replace wiht a dummy year
  mutate(year = str_replace(season,"(.+)/(.+)","\\1") %>% as.double())

denv %>% glimpse()

## Rows: 468
## Columns: 12
## $ season               <chr> "2000/2001", "2000/2001", "2000/2001", "2000/2001~
## $ season_week          <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15~
## $ week_start_date      <date> 2000-07-01, 2000-07-08, 2000-07-15, 2000-07-22, ~
## $ denv1_cases          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
## $ denv2_cases          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
## $ denv3_cases          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
## $ denv4_cases          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
## $ other_positive_cases <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1~
## $ total_cases          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1~
## $ year                 <dbl> 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2~
## $ epiweek              <dbl> 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 3~
## $ adm                  <chr> "iquitos", "iquitos", "iquitos", "iquitos", "iqui~

datos de población

popdb <-
  readr::read_csv("https://predict.cdc.gov/api/v1/attachments/dengue%20forecasting%20project/iquitos_population_data.csv") %>%
  janitor::clean_names() %>%
  mutate(adm="iquitos")

popdb %>% glimpse()

## Rows: 15
## Columns: 3
## $ year                 <dbl> 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2~
## $ estimated_population <dbl> 386666, 393355, 399770, 405988, 412095, 418168, 4~
## $ adm                  <chr> "iquitos", "iquitos", "iquitos", "iquitos", "iqui~

adaptar bases de datos

Segundo, adaptar ambas bases de datos con la función epi_adapt_timeserie()

epi_adapted <-
  epi_adapt_timeserie(db_disease = denv,
                      db_population = popdb,
                      var_admx = adm,
                      var_year = year, # must be a common variable name between datasets
                      var_week = season_week,
                      # var_year = year, 
                      # var_week = epiweek,
                      var_event_count = total_cases,
                      var_population = estimated_population)

## Joining, by = c("var_admx", "var_year")

# epi_adapted

filtrar por años

Tercero, filtrar por años para discriminar entre datos de años históricos (diferentes al último año) y año corriente (iguales al último año) con la función filter().
Filtrar filas usando el paquete dplyr (tutorial primer)

disease_now <- epi_adapted %>%
  filter(var_year==max(var_year))

disease_pre <- epi_adapted %>%
  filter(var_year!=max(var_year))

crear canal endémico

Cuarto, crear el canal endémico con la función epi_create_channel()
- Aquí puedes elegir entre tres métodos disponibles (Bortman, 1999):
  - "gmean_1sd" es media geométrica con 1 desviación estándar (por defecto).
  - "gmean_2sd" es media geométrica con 2 desviaciones estándar.
  - "gmean_ci" es media geométrica con intervalos de confianza al 95%.

disease_channel <-
  epi_create_channel(time_serie = disease_pre,
                     disease_name = "denv",
                     method = "gmean_1sd")

## Joining, by = "var_admx"

disease_channel

## # A tibble: 52 x 6
##    var_admx var_week median  low_95 upp_95 key  
##    <fct>       <dbl>  <dbl>   <dbl>  <dbl> <chr>
##  1 iquitos         1   2.80  0.267    6.80 denv 
##  2 iquitos         2   2.64  0.484    5.82 denv 
##  3 iquitos         3   2.78  0.259    6.75 denv 
##  4 iquitos         4   2.86  0.424    6.62 denv 
##  5 iquitos         5   2.72  0.454    6.12 denv 
##  6 iquitos         6   2.25  0.127    5.44 denv 
##  7 iquitos         7   2.22  0.0474   5.52 denv 
##  8 iquitos         8   2.34  0.167    5.61 denv 
##  9 iquitos         9   1.82 -0.420    5.44 denv 
## 10 iquitos        10   2.66 -0.513    8.65 denv 
## # ... with 42 more rows

unir y graficar canal

Finalmente, unir bases con epi_join_channel() y
Graficar usando ggplot con epi_plot_channel()

epi_join_channel(disease_channel = disease_channel,
                 disease_now = disease_now) %>%
  # ggplot
  epi_plot_channel() +
  labs(title = "Dengue virus Endemic Channel. Iquitos, Peru 2008/2009",
       caption = "Source: https://dengueforecasting.noaa.gov/",
       # x = "epiweeks",
       x = "Seasonal week",
       y = "Number of cases") +
  theme_bw()

## Joining, by = c("var_admx", "var_week")

P4.1: ¿Por qué crees que hay una disminución repentina de casos y de zonas del canal entre las semanas estacionales 20 y 30 (diciembre - enero)?

P4.2: ¿Preguntas?

CODA

05 - difusión

Aquí lista de recursos útiles para construir y publicar tableros de comando o dashboards

Visualizando datos de Salud Pública y Epidemiología de campo para el Diplomado 2021 en diapositivas
LatinR 2020: Aplicaciones web interactivas con Shiny - video y repositorio de materiales
Rstudio Webinar sobre el paquete flexdashboards

06 - cierre

¡Gracias por su atención!

Andree Valle Campos

`@avallecam` click aquí

avallecam@gmail.com

código completo

# instalar paquetes si se requiere ----------------------------------------

if(!require("remotes")) install.packages("remotes")
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("lubridate")) install.packages("lubridate")
if(!require("compareGroups")) install.packages("compareGroups")
if(!require("rio")) install.packages("rio")
if(!require("naniar")) install.packages("naniar")
if(!require("incidence")) install.packages("incidence")
if(!require("epiR")) install.packages("epiR")
if(!require("epiDisplay")) install.packages("epiDisplay")
if(!require("epitidy")) remotes::install_github("avallecam/epitidy")
if(!require("incidenceflow")) remotes::install_github("avallecam/incidenceflow")
if(!require("epichannel")) remotes::install_github("avallecam/epichannel")


# # 01 - tiempo espacio persona ------------------------------------
# ## _librerias a usar ---------------------------------

library(tidyverse)
library(lubridate)
library(compareGroups)

# ## _importar base de datos ---------------------------------

ruta_limpio <- "https://github.com/avallecam/epiapli2019/raw/master/data/casoslimpio_20190916.rds"
casos_limpio <- rio::import(file = ruta_limpio)
casos_limpio

casos_limpio %>% naniar::miss_var_summary()
casos_limpio %>% naniar::vis_miss()

# ## _distribución en tiempo ---------------------------------

casos_limpio %>%
  ggplot(aes(x = date_of_onset)) +
  geom_histogram(binwidth = 7, color="white") +
  scale_x_date(date_breaks = "7 day", date_labels = "%b-%d") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# ## _distribución en espacio ---------------------------------

casos_limpio %>%
  ggplot(aes(x = lon, y = lat, colour = date_of_onset)) +
  geom_point() +
  scale_color_gradient(low = "red", high = "yellow", trans = "date") +
  theme_bw()

# ## _tabla descriptiva ---------------------------------

casos_limpio %>%
  select(outcome, gender, hospital) %>%
  mutate(hospital = fct_infreq(hospital),
         outcome = fct_infreq(outcome)) %>%
  compareGroups(~., data = .) %>%
  createTable() %>%
  export2md()


# # 02 - factores asociados ------------------------------------
# ## _librerias a usar ---------------------------------

library(epitidy)

# ## _importar base de datos ---------------------------------

## # dentro de R
## if(!require("mosaicData")) install.packages("mosaicData")
## ?mosaicData::Whickham

smoke_limpio <- "https://github.com/avallecam/epiapli2019/raw/master/data/smokeclean_20190906.rds"
smoke_clean <- rio::import(file = smoke_limpio)
smoke_clean

# ## _tabla descriptiva ---------------------------------

tabla1 <- smoke_clean %>%
  compareGroups(outcome~smoker+agegrp+age, data = ., byrow = T) %>%
  createTable(show.all = T, sd.type = 2)

tabla1

## tabla1 %>%
##   export2xls("tabla/tabla02.xls")

# ## _medidas de asociación ---------------------------------

# ###__ modelo simple ---------------------------------

# smoke_clean %>% glimpse()

#simple
simple_model <- glm(outcome_1 ~ smoker,
                    data = smoke_clean,
                    family = poisson(link = "log"))
epi_tidymodel_rr(simple_model)



# ###__ modelo múltiple ---------------------------------

#multiple: controlar por confusión
multiple_model <- glm(outcome_1 ~ smoker + age,
                      data = smoke_clean,
                      family = poisson(link = "log"))
epi_tidymodel_rr(multiple_model)

## epi_tidymodel_rr(multiple_model) %>%
##   writexl::write_xlsx("tabla/tabla03.xlsx")

## # medidas de asociación ---------------------------------------------------
##
## #epiR
## library(epiR)
## smoke_tabla1 <- with(smoke_clean,table(smoker_2,outcome_2)) %>% print()
## epi.2by2(smoke_tabla1,method = "cohort.count")
##
## #epiDisplay
## library(epiDisplay)
## smoke_tab2 <- with(smoke_clean,table(outcome,smoker)) %>% print()
## cs(cctable = smoke_tab2)
##
## # controlar por confusión -------------------------------------------------
##
## #Mantel-Haenszel
## smoke_tab3 <- with(smoke_clean,table(smoker_2,outcome_2,agegrp)) %>% print()
## epi.2by2(smoke_tab3,method = "cohort.count")
## mhor(mhtable=smoke_tab3,graph = F,design = "cohort")


# # 03 - curva epidémica ------------------------------------
# ## _librerias a usar ---------------------------------

library(outbreaks) #sample data
library(incidence) #core functions
library(incidenceflow)

# ## _importar base de datos ---------------------------------

ebola_sim$linelist %>% as_tibble()

dat <- ebola_sim$linelist$date_of_onset
enframe(dat,value = "date_of_onset") %>% select(date_of_onset)

i.7 <- incidence(dat, interval=7)
i.7

plot(i.7)

# ## _curva ascendente ---------------------------------

f1 <- fit(i.7[1:20])
f1

# ## _curva descendente ---------------------------------

f2 <- fit_optim_split(i.7)
f2

# ## _generar tablas resumen ---------------------------------
# ###__ curva ascendente ---------------------------------

f1 %>% tidy_incidence()

f1 %>% glance_incidence()

# ###__ curva descendente ---------------------------------

f2 %>% pluck("fit") %>% tidy_incidence()

f2 %>% pluck("fit") %>% glance_incidence()

# # 04 - canal endémico ------------------------------------
# ## _librerias a usar ---------------------------------

library(epichannel)

# ## _importar base de datos ---------------------------------
# ###__ datos de enfermedad ---------------------------------

denv <-
  readr::read_csv("https://predict.cdc.gov/api/v1/attachments/dengue%20forecasting%20project/iquitos_training_data.csv") %>%
  mutate(year = lubridate::year(week_start_date),
         epiweek = lubridate::epiweek(week_start_date)) %>%
  mutate(adm="iquitos") %>%
  # cases per season - replace wiht a dummy year
  mutate(year = str_replace(season,"(.+)/(.+)","\\1") %>% as.double())

denv %>% glimpse()

# ###__ datos de población ---------------------------------

popdb <-
  readr::read_csv("https://predict.cdc.gov/api/v1/attachments/dengue%20forecasting%20project/iquitos_population_data.csv") %>%
  janitor::clean_names() %>%
  mutate(adm="iquitos")

popdb %>% glimpse()

## # denv %>% count(year,season,lag_year)
## # denv %>%
## #   ggplot(aes(x = week_start_date,y = total_cases)) +
## #   geom_col()
## # popdb %>% count(year)
## # denv %>% count(year)
## # denv %>% left_join(popdb)

# ## _adaptar bases de datos ---------------------------------

epi_adapted <-
  epi_adapt_timeserie(db_disease = denv,
                      db_population = popdb,
                      var_admx = adm,
                      var_year = year, # must be a common variable name between datasets
                      var_week = season_week,
                      # var_year = year,
                      # var_week = epiweek,
                      var_event_count = total_cases,
                      var_population = estimated_population)
epi_adapted

# ## _filtrar por años ---------------------------------

disease_now <- epi_adapted %>%
  filter(var_year==max(var_year))

disease_pre <- epi_adapted %>%
  filter(var_year!=max(var_year))

# ## _crear canal endémico ---------------------------------

disease_channel <-
  epi_create_channel(time_serie = disease_pre,
                     disease_name = "denv",
                     method = "gmean_1sd")
disease_channel

# ## _unir y graficar canal ---------------------------------

epi_join_channel(disease_channel = disease_channel,
                 disease_now = disease_now) %>%
  # ggplot
  epi_plot_channel() +
  labs(title = "Dengue virus Endemic Channel. Iquitos, Peru 2008/2009",
       caption = "Source: https://dengueforecasting.noaa.gov/",
       # x = "epiweeks",
       x = "Seasonal week",
       y = "Number of cases") +
  theme_bw()

Vigilancia en Salud Pública: Análisis de datos y difusión de la información

Andree Valle Campos avallecam@gmail.com

2021-05-31

Advertencias

Objetivos

INICIO

00 - ¿cómo empezar con R?

paso a paso

algunos conceptos clave

01 - tiempo espacio persona

librerias a usar

importar base de datos

distribución en tiempo

distribución en espacio

tabla descriptiva

02 - factores asociados

librerias a usar

importar base de datos

tabla descriptiva

medidas de asociación

modelo simple

modelo múltiple

métodos alternativo

INTERMEDIO

03 - curva epidémica

librerias a usar

importar base de datos

curva ascendente

curva descendente

generar tablas resumen

curva ascendente

curva descendente

04 - canal endémico

librerias a usar

importar base de datos

datos de enfermedad

datos de población

adaptar bases de datos

filtrar por años

crear canal endémico

unir y graficar canal

CODA

05 - difusión

06 - cierre

¡Gracias por su atención!

Andree Valle Campos @avallecam click aquí avallecam@gmail.com

código completo

Andree Valle Campos
avallecam@gmail.com

Andree Valle Campos

`@avallecam` click aquí

avallecam@gmail.com