Análisis de brotes en tiempo real: el ébola como estudio de caso
Anne Cori, Natsuko Imai, Finlay Campbell, Zhian N. Kamvar, Thibaut Jombart,José M. Velasco-España, Cándida Díaz-Brochero, Zulma M. Cucunubá
2022-10-25
Introducción
Un nuevo brote de virus del Ébola (EVE) en un país ficticio de África occidental
Objetivos:
-Identificar los parámetros necesarios en casos de transmisión de enfermedades infecciosas de persona a persona. -Estimar e interpretar la tasa de crecimiento y el tiempo en que se duplica la epidemia. -Estimar el intervalo serial a partir de los datos pareados de individuos infectantes/ individuos infectados. -Estimar e interpretar el número de reproducción instantáneo de la epidemia -Estimar la probabilidad de letalidad (CFR) -Calcular y graficar la incidencia
Conceptos básicos a desarrollar
En esta práctica se desarrollarán los siguientes conceptos:
- Transmisión de enfermedades infecciosas de persona a persona
- Número reproductivo efectivo
- Probabilidad de letalidad
- Intervalo serial
- Tasa de crecimiento
- Incidencia
Cargue de librerías:
library(readxl)
library(incidence)
library(epicontacts)
library(distcrete)
library(epitrix)
library(EpiEstim)
library(projections)
library(ggplot2)
library(magrittr)
library(binom)
library(ape)
library(outbreaker2)
library(tidyverse)
library(knitr)Cargue de bases de datos
linelist_clean <- readRDS("data/linelist_clean.rds")
contacts <- read_excel("data/contacts_20140701.xlsx", na = c("", "NA"))Estructura de los datos
str(contacts)tibble [60 × 3] (S3: tbl_df/tbl/data.frame)
$ infector: chr [1:60] "d1fafd" "f5c3d8" "0f58c4" "f5c3d8" ...
$ case_id : chr [1:60] "53371b" "0f58c4" "881bd4" "d58402" ...
$ source : chr [1:60] "other" "other" "other" "other" ...str(linelist_clean)tibble [166 × 11] (S3: tbl_df/tbl/data.frame)
$ case_id : chr [1:166] "d1fafd" "53371b" "f5c3d8" "6c286a" ...
$ generation : num [1:166] 0 1 1 2 2 0 3 3 2 3 ...
$ date_of_infection : Date[1:166], format: NA "2014-04-09" ...
$ date_of_onset : Date[1:166], format: "2014-04-07" "2014-04-15" ...
$ date_of_hospitalisation: Date[1:166], format: "2014-04-17" "2014-04-20" ...
$ date_of_outcome : Date[1:166], format: "2014-04-19" NA ...
$ outcome : Factor w/ 2 levels "Death","Recover": NA NA 2 1 2 NA 2 1 2 1 ...
$ gender : Factor w/ 2 levels "f","m": 1 2 1 1 1 1 1 1 2 2 ...
$ hospital : Factor w/ 6 levels "Connaught Hospital",..: 2 1 4 NA 4 NA 1 5 4 6 ...
$ lon : num [1:166] -13.2 -13.2 -13.2 -13.2 -13.2 ...
$ lat : num [1:166] 8.47 8.46 8.48 8.46 8.45 ...#Mire más de cerca los datos contenidos en este
linelist.
head(linelist_clean)# A tibble: 6 × 11
case_id generation date_of_i…¹ date_of_…² date_of_…³ date_of_…⁴ outcome gender
<chr> <dbl> <date> <date> <date> <date> <fct> <fct>
1 d1fafd 0 NA 2014-04-07 2014-04-17 2014-04-19 <NA> f
2 53371b 1 2014-04-09 2014-04-15 2014-04-20 NA <NA> m
3 f5c3d8 1 2014-04-18 2014-04-21 2014-04-25 2014-04-30 Recover f
4 6c286a 2 NA 2014-04-27 2014-04-27 2014-05-07 Death f
5 0f58c4 2 2014-04-22 2014-04-26 2014-04-29 2014-05-17 Recover f
6 49731d 0 2014-03-19 2014-04-25 2014-05-02 2014-05-07 <NA> f
# … with 3 more variables: hospital <fct>, lon <dbl>, lat <dbl>, and
# abbreviated variable names ¹date_of_infection, ²date_of_onset,
# ³date_of_hospitalisation, ⁴date_of_outcomeProbabilidad de muerte en los casos reportados (CFR)
table(linelist_clean$outcome, useNA = "ifany")
Death Recover <NA>
60 43 63 n_dead <- sum(linelist_clean$outcome %in% "Death")
n_known_outcome <- sum(linelist_clean$outcome %in% c("Death", "Recover"))
cfr <- n_dead / n_known_outcome
cfr_with_CI <- binom.confint(n_dead, n_known_outcome, method = "exact")
kable(cfr_with_CI, caption = "cfr_with_CI")| method | x | n | mean | lower | upper |
|---|---|---|---|---|---|
| exact | 60 | 103 | 0.5825243 | 0.4812264 | 0.6789504 |
Curvas de incidencia
i_daily <- incidence(linelist_clean$date_of_onset)
i_daily<incidence object>
[166 cases from days 2014-04-07 to 2014-06-29]
$counts: matrix with 84 rows and 1 columns
$n: 166 cases in total
$dates: 84 dates marking the left-side of bins
$interval: 1 day
$timespan: 84 days
$cumulative: FALSEplot(i_daily, border = "black")
Cálculo de la incidencia semanal i_weekly
i_weekly <- incidence(linelist_clean$date_of_onset, interval = 7,
last_date = as.Date(max(linelist_clean$date_of_hospitalisation, na.rm = TRUE)))
i_weekly<incidence object>
[166 cases from days 2014-04-07 to 2014-06-30]
[166 cases from ISO weeks 2014-W15 to 2014-W27]
$counts: matrix with 13 rows and 1 columns
$n: 166 cases in total
$dates: 13 dates marking the left-side of bins
$interval: 7 days
$timespan: 85 days
$cumulative: FALSEplot(i_weekly, border = "black")
Estimación de la tasa de crecimiento mediante un modelo log-lineal
Grafique la incidencia transformada logarítmicamente:
ggplot(as.data.frame(i_weekly)) +
geom_point(aes(x = dates, y = log(counts))) +
scale_x_incidence(i_weekly) +
xlab("date") +
ylab("log weekly incidence") +
theme_minimal()
#### Ajuste un modelo log-lineal a los datos de incidencia semanal
f <- incidence::fit(i_weekly)
f<incidence_fit object>
$model: regression of log-incidence over time
$info: list containing the following items:
$r (daily growth rate):
[1] 0.04145251
$r.conf (confidence interval):
2.5 % 97.5 %
[1,] 0.02582225 0.05708276
$doubling (doubling time in days):
[1] 16.72148
$doubling.conf (confidence interval):
2.5 % 97.5 %
[1,] 12.14285 26.84302
$pred: data.frame of incidence predictions (12 rows, 5 columns)plot(i_weekly, fit = f)
Encontrando una fecha límite adecuada para el modelo log-lineal, en función de los retrasos observados
Es posible que desee examinar cuánto tiempo después de la aparición de los síntomas los casos son hospitalizados
summary(as.numeric(linelist_clean$date_of_hospitalisation - linelist_clean$date_of_onset)) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 1.00 2.00 3.53 5.00 22.00 Semanas a descartar al final de la epicurva
n_weeks_to_discard <- 2
min_date <- min(i_daily$dates)
max_date <- max(i_daily$dates) - n_weeks_to_discard * 7
# Para truncar la incidencia semanal
i_weekly_trunc <- subset(i_weekly,
from = min_date,
to = max_date) # descarte las últimas semanas de datos
# incidencia diaria truncada.No la usamos para la regresión lineal pero se puede usar más adelante
i_daily_trunc <- subset(i_daily,
from = min_date,
to = max_date) # eliminamos las últimas dos semanas de datosVuelva a montar y a graficar el modelo logarítmico lineal, pero
utilizando los datos truncados i_weekly_trunc.
f <- incidence::fit(i_weekly_trunc)
f<incidence_fit object>
$model: regression of log-incidence over time
$info: list containing the following items:
$r (daily growth rate):
[1] 0.05224047
$r.conf (confidence interval):
2.5 % 97.5 %
[1,] 0.03323024 0.0712507
$doubling (doubling time in days):
[1] 13.2684
$doubling.conf (confidence interval):
2.5 % 97.5 %
[1,] 9.728286 20.85893
$pred: data.frame of incidence predictions (10 rows, 5 columns)plot(i_weekly_trunc, fit = f)
Observe las estadísticas resumidas de su ajuste:
summary(f$model)
Call:
stats::lm(formula = log(counts) ~ dates.x, data = df)
Residuals:
Min 1Q Median 3Q Max
-0.73474 -0.31655 -0.03211 0.41798 0.65311
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.186219 0.332752 0.560 0.591049
dates.x 0.052240 0.008244 6.337 0.000224 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5241 on 8 degrees of freedom
Multiple R-squared: 0.8339, Adjusted R-squared: 0.8131
F-statistic: 40.16 on 1 and 8 DF, p-value: 0.0002237Estimacion de la tasa de crecimiento
daily_growth_rate <- f$model$coefficients['dates.x']
daily_growth_rate dates.x
0.05224047 # intervalo de confianza:
daily_growth_rate_CI <- confint(f$model, 'dates.x', level=0.95)
daily_growth_rate_CI 2.5 % 97.5 %
dates.x 0.03323024 0.0712507Estimacion del tiempo de duplicacion
doubling_time_days <- log(2) / daily_growth_rate
doubling_time_days dates.x
13.2684 # intervalo de confianza:
doubling_time_days_CI <- log(2) / rev(daily_growth_rate_CI)
doubling_time_days_CI[1] 9.728286 20.858930Rastreo de contactos
Generacion de la red de rastreo de contactos:
Usando la función make_epicontacts en el paquete
epicontacts, cree un nuevo objeto llamado
epi_contacts . Asegúrese de comprobar los nombres de las
columnas de los argumentos relevantes “to” y
“from”.
epi_contacts <- make_epicontacts(linelist_clean,
contacts,
id = "case_id",
from = "infector",
to = "case_id")
epi_contacts
/// Epidemiological Contacts //
// class: epicontacts
// 166 cases in linelist; 60 contacts; non directed
// linelist
# A tibble: 166 × 11
id generation date_of_i…¹ date_of_…² date_of_…³ date_of_…⁴ outcome gender
<chr> <dbl> <date> <date> <date> <date> <fct> <fct>
1 d1fafd 0 NA 2014-04-07 2014-04-17 2014-04-19 <NA> f
2 53371b 1 2014-04-09 2014-04-15 2014-04-20 NA <NA> m
3 f5c3d8 1 2014-04-18 2014-04-21 2014-04-25 2014-04-30 Recover f
4 6c286a 2 NA 2014-04-27 2014-04-27 2014-05-07 Death f
5 0f58c4 2 2014-04-22 2014-04-26 2014-04-29 2014-05-17 Recover f
6 49731d 0 2014-03-19 2014-04-25 2014-05-02 2014-05-07 <NA> f
7 f9149b 3 NA 2014-05-03 2014-05-04 NA Recover f
8 881bd4 3 2014-04-26 2014-05-01 2014-05-05 NA Death f
9 e66fa4 2 NA 2014-04-21 2014-05-06 2014-05-12 Recover m
10 20b688 3 NA 2014-05-05 2014-05-06 2014-05-14 Death m
# … with 156 more rows, 3 more variables: hospital <fct>, lon <dbl>, lat <dbl>,
# and abbreviated variable names ¹date_of_infection, ²date_of_onset,
# ³date_of_hospitalisation, ⁴date_of_outcome
// contacts
# A tibble: 60 × 3
from to source
<chr> <chr> <chr>
1 d1fafd 53371b other
2 f5c3d8 0f58c4 other
3 0f58c4 881bd4 other
4 f5c3d8 d58402 other
5 20b688 d8a13d funeral
6 2ae019 a3c8b8 other
7 20b688 974bc1 funeral
8 2ae019 72b905 funeral
9 40ae5f b8f2fd funeral
10 f1f60f 09e386 other
# … with 50 more rows# observe la fuente de infección reportada de los contactos.
table(epi_contacts$contacts$source, useNA = "ifany")##
## funeral other
## 20 40p <- plot(epi_contacts, node_shape = "gender", shapes = c(m = "male", f = "female"), node_color = "gender", edge_color = "source", selector = FALSE)
pUsando la función match ( ver? Match )
verifique que los contactos visualizados sean realmente casos.
match(contacts$case_id, linelist_clean$case_id)## [1] 2 5 8 14 15 16 18 20 21 22 24 25 26 27 30 33 37 38 40
## [20] NA 47 50 57 58 61 62 66 67 68 69 71 73 77 82 84 86 88 92
## [39] 93 94 96 101 106 112 113 119 123 128 130 139 142 143 144 145 150 152 154
## [58] 157 159 163Estimación del intervalo serial (SI)
si_obs <- get_pairwise(epi_contacts, "date_of_onset")
summary(si_obs) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
2.000 5.000 7.000 9.254 12.500 30.000 1 hist(si_obs, breaks = 0:30,
xlab = "Días después de la aparición de los síntomas", ylab = "Frecuencia",
main = "intervalo serial (distribución empírica)",
col = "grey", border = "white")
Ajuste a distribucion gamma
si_fit <- fit_disc_gamma(si_obs, w = 1)
si_fit$mu
[1] 8.752409
$cv
[1] 0.6963483
$sd
[1] 6.094725
$ll
[1] -179.6802
$converged
[1] TRUE
$distribution
A discrete distribution
name: gamma
parameters:
shape: 2.0622770570919
scale: 4.24405104660505#comparacion de la distribucion ajustada
hist(si_obs, xlab = "Días después de la aparición de los síntomas", ylab = "Frecuencia",
main = "intervalo serial: ajustar a los datos", col = "salmon", border = "white",
50, ylim = c(0, 0.15), freq = FALSE, breaks = 0:35)
# Por favor revisar aquí como sacar el fit porque no está funcionando. Gracias, Zulma.
# points(0:60, si_fit$d(0:60), col = "#9933ff", pch = 20)
# points(0:60, si_fit$sd(0:60), col = "#9933ff", type = "l", lty = 2)Estimación de la transmisibilidad variable en el tiempo, R(t)
Cuando la suposición de que (\(R\))
es constante en el tiempo se vuelve insostenible, una alternativa es
estimar la transmisibilidad variable en el tiempo utilizando el número
de reproducción instantánea (\(R_t\)).
Este enfoque, introducido por Cori et al. (2013), se implementa en el
paquete EpiEstim. Estima (\(R_t\)) para ventanas de tiempo
personalizadas (el valor predeterminado es una sucesión de ventanas de
tiempo deslizantes), utilizando la probabilidad de Poisson. A
continuación, estimamos la transmisibilidad para ventanas de tiempo
deslizantes de 1 semana (el valor predeterminado de
estimate_R):
config <- make_config(mean_si = si_fit$mu, # media de la distribución si estimada anteriormente
std_si = si_fit$sd, # desviación estándar de la distribución si estimada anteriormente
t_start = 2, # día de inicio de la ventana de tiempo
t_end = length(i_daily_trunc$counts)) # último día de la ventana de tiempoconfig = make_config(list(mean_si = si_fit$mu, std_si = si_fit$sd))
# t_start y t_end se configuran automáticamente para estimar R en ventanas deslizantes para 1 semana de forma predeterminada.# use estimate_R using method = "parametric_si"
Rt <- estimate_R(i_daily_trunc, method = "parametric_si",
si_data = si_data,
config = config)
# mire las estimaciones de Rt más recientes:
tail(Rt$R[, c("t_start", "t_end", "Median(R)",
"Quantile.0.025(R)", "Quantile.0.975(R)")]) t_start t_end Median(R) Quantile.0.025(R) Quantile.0.975(R)
58 59 65 1.3243809 0.8686233 1.917253
59 60 66 1.4470115 0.9730882 2.054542
60 61 67 1.2504921 0.8201617 1.810287
61 62 68 1.0120489 0.6365492 1.512537
62 63 69 0.8445741 0.5099601 1.301134
63 64 70 1.0282900 0.6543105 1.523422Grafique la estimación de \(R\) sobre le tiempo:
plot(Rt, legend = FALSE)
Guardar salidas
saveRDS(epi_contacts, "data/epi_contacts.rds")
saveRDS(si_fit, "data/si_fit.rds")Sobre este documento
Este documento ha sido una adaptación de los materiales originales disponibles en RECON Learn
Contribuciones
Autores originales:
Anne Cori
Natsuko Imai
Finlay Campbell
Zhian N. Kamvar
Thibaut Jombart
Cambios menores y adaptación a español:
José M. Velasco-España
Cándida Díaz-Brochero
Zulma M. Cucunubá