-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathEnfermedadX2.Rmd
1337 lines (985 loc) · 51.9 KB
/
EnfermedadX2.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
title: "Taller Día 4 - Grupo 2 - Estimación de las distribuciones de rezagos epidemiológicos: Enfermedad X"
author: "Kelly Charniga, PhD MPH & Zulma Cucunubá MD, PhD"
date: "2023-11-23"
output: html_document
teaching: 90
exercises: 8
---
```{r include=FALSE}
## Calcule la verosimilitud DIC mediante integración
diclik <- function(par1, par2, EL, ER, SL, SR, dist){
## Si la ventana de síntomas es mayor que la ventana de exposición
if(SR-SL>ER-EL){
dic1 <- integrate(fw1, lower=SL-ER, upper=SL-EL,
subdivisions=10,
par1=par1, par2=par2,
EL=EL, ER=ER, SL=SL, SR=SR,
dist=dist)$value
if (dist == "W"){
dic2 <- (ER-EL)*
(pweibull(SR-ER, shape=par1, scale=par2) - pweibull(SL-EL, shape=par1, scale=par2))
} else if (dist == "off1W"){
dic2 <- (ER-EL)*
(pweibullOff1(SR-ER, shape=par1, scale=par2) - pweibullOff1(SL-EL, shape=par1, scale=par2))
} else if (dist == "G"){
dic2 <- (ER-EL)*
(pgamma(SR-ER, shape=par1, scale=par2) - pgamma(SL-EL, shape=par1, scale=par2))
} else if (dist == "off1G"){
dic2 <- (ER-EL)*
(pgammaOff1(SR-ER, shape=par1, scale=par2) - pgammaOff1(SL-EL, shape=par1, scale=par2))
} else if (dist == "L") {
dic2 <- (ER-EL)*
(plnorm(SR-ER, par1, par2) - plnorm(SL-EL, par1, par2))
} else if (dist == "off1L") {
dic2 <- (ER-EL)*
(plnormOff1(SR-ER, par1, par2) - plnormOff1(SL-EL, par1, par2))
} else {
stop("distribution not supported")
}
dic3 <- integrate(fw3, lower=SR-ER, upper=SR-EL,
subdivisions=10,
par1=par1, par2=par2,
EL=EL, ER=ER, SL=SL, SR=SR,
dist=dist)$value
return(dic1 + dic2 + dic3)
}
## Si la ventana de exposición es mayor que la ventana de síntomas
else{
dic1 <- integrate(fw1, lower=SL-ER, upper=SR-ER, subdivisions=10,
par1=par1, par2=par2,
EL=EL, ER=ER, SL=SL, SR=SR,
dist=dist)$value
if (dist == "W"){
dic2 <- (SR-SL)*
(pweibull(SL-EL, shape=par1, scale=par2) - pweibull(SR-ER, shape=par1, scale=par2))
} else if (dist == "off1W"){
dic2 <- (SR-SL)*
(pweibullOff1(SL-EL, shape=par1, scale=par2) - pweibullOff1(SR-ER, shape=par1, scale=par2))
} else if (dist == "G"){
dic2 <- (SR-SL)*
(pgamma(SL-EL, shape=par1, scale=par2) - pgamma(SR-ER, shape=par1, scale=par2))
} else if (dist == "off1G"){
dic2 <- (SR-SL)*
(pgammaOff1(SL-EL, shape=par1, scale=par2) - pgammaOff1(SR-ER, shape=par1, scale=par2))
} else if (dist == "L"){
dic2 <- (SR-SL)*
(plnorm(SL-EL, par1, par2) - plnorm(SR-ER, par1, par2))
} else if (dist == "off1L"){
dic2 <- (SR-SL)*
(plnormOff1(SL-EL, par1, par2) - plnormOff1(SR-ER, par1, par2))
} else {
stop("distribution not supported")
}
dic3 <- integrate(fw3, lower=SL-EL, upper=SR-EL,
subdivisions=10,
par1=par1, par2=par2,
EL=EL, ER=ER, SL=SL, SR=SR,
dist=dist)$value
return(dic1 + dic2 + dic3)
}
}
## Esta verosimilitud DIC está diseñada para datos que tienen intervalos superpuestos
diclik2 <- function(par1, par2, EL, ER, SL, SR, dist){
if(SL>ER) {
return(diclik(par1, par2, EL, ER, SL, SR, dist))
} else {
lik1 <- integrate(diclik2.helper1, lower=EL, upper=SL,
SL=SL, SR=SR, par1=par1, par2=par2, dist=dist)$value
lik2 <- integrate(diclik2.helper2, lower=SL, upper=ER,
SR=SR, par1=par1, par2=par2, dist=dist)$value
return(lik1+lik2)
}
}
## Funciones de verosimilitud para diclik2
diclik2.helper1 <- function(x, SL, SR, par1, par2, dist){
if (dist =="W"){
pweibull(SR-x, shape=par1, scale=par2) - pweibull(SL-x, shape=par1, scale=par2)
} else if (dist =="off1W") {
pweibullOff1(SR-x, shape=par1, scale=par2) - pweibullOff1(SL-x, shape=par1, scale=par2)
} else if (dist =="G") {
pgamma(SR-x, shape=par1, scale=par2) - pgamma(SL-x, shape=par1, scale=par2)
} else if (dist=="off1G"){
pgammaOff1(SR-x, shape=par1, scale=par2) - pgammaOff1(SL-x, shape=par1, scale=par2)
} else if (dist == "L"){
plnorm(SR-x, par1, par2) - plnorm(SL-x, par1, par2)
} else if (dist == "off1L"){
plnormOff1(SR-x, par1, par2) - plnormOff1(SL-x, par1, par2)
} else {
stop("distribution not supported")
}
}
diclik2.helper2 <- function(x, SR, par1, par2, dist){
if (dist =="W"){
pweibull(SR-x, shape=par1, scale=par2)
} else if (dist =="off1W") {
pweibullOff1(SR-x, shape=par1, scale=par2)
} else if (dist =="G") {
pgamma(SR-x, shape=par1, scale=par2)
} else if (dist =="off1G") {
pgammaOff1(SR-x, shape=par1, scale=par2)
} else if (dist=="L"){
plnorm(SR-x, par1, par2)
} else if (dist=="off1L"){
plnormOff1(SR-x, par1, par2)
} else {
stop("distribution not supported")
}
}
## Funciones que manipulan/calculan la verosimilitud para los datos censurados
## Las funciones codificadas aquí se toman directamente de las
## notas de verosimilitud censurada por intervalos dobles.
fw1 <- function(t, EL, ER, SL, SR, par1, par2, dist){
## Función que calcula la primera función para la integral DIC
if (dist=="W"){
(ER-SL+t) * dweibull(x=t,shape=par1,scale=par2)
} else if (dist=="off1W") {
(ER-SL+t) * dweibullOff1(x=t,shape=par1,scale=par2)
} else if (dist=="G") {
(ER-SL+t) * dgamma(x=t, shape=par1, scale=par2)
} else if (dist=="off1G") {
(ER-SL+t) * dgammaOff1(x=t, shape=par1, scale=par2)
} else if (dist =="L"){
(ER-SL+t) * dlnorm(x=t, meanlog=par1, sdlog=par2)
} else if (dist =="off1L"){
(ER-SL+t) * dlnormOff1(x=t, meanlog=par1, sdlog=par2)
} else {
stop("distribution not supported")
}
}
fw3 <- function(t, EL, ER, SL, SR, par1, par2, dist){
## Función que calcula la tercera función para la integral DIC
if (dist == "W"){
(SR-EL-t) * dweibull(x=t, shape=par1, scale=par2)
} else if (dist == "off1W"){
(SR-EL-t) * dweibullOff1(x=t, shape=par1, scale=par2)
} else if (dist == "G"){
(SR-EL-t) * dgamma(x=t, shape=par1, scale=par2)
} else if (dist == "off1G"){
(SR-EL-t) * dgammaOff1(x=t, shape=par1, scale=par2)
} else if (dist == "L") {
(SR-EL-t) * dlnorm(x=t, meanlog=par1, sdlog=par2)
} else if (dist == "off1L"){
(SR-EL-t) * dlnormOff1(x=t, meanlog=par1, sdlog=par2)
} else {
stop("distribution not supported")
}
}
```
:::::::::::::::::::::::::::::::::::::: questions
- ¿Cómo responder ante un brote de una enfermedad desconocida?
::::::::::::::::::::::::::::::::::::::::::::::::
::::::::::::::::::::::::::::::::::::: objectives
Al final de este taller usted podrá:
- Comprender los conceptos clave de las distribuciones de rezagos epidemiológicos para la Enfermedad X.
- Comprender las estructuras de datos y las herramientas para el análisis de datos de rastreo de contactos.
- Aprender a ajustar las estimaciones del intervalo serial y el período de incubación de la Enfermedad X teniendo en cuenta la censura por intervalo utilizando un marco de trabajo Bayesiano.
- Aprender a utilizar estos parámetros para informar estrategias de control en un brote de un patógeno desconocido.
::::::::::::::::::::::::::::::::::::::::::::::::
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## 1. Introducción
La Enfermedad X representa un hipotético, pero plausible, brote de una enfermedad infecciosa en el futuro. Este término fue acuñado por la Organización Mundial de la Salud (OMS) y sirve como un término general para un patógeno desconocido que podría causar una epidemia grave a nivel internacional. Este concepto representa la naturaleza impredecible de la aparición de enfermedades infecciosas y resalta la necesidad de una preparación global y mecanismos de respuesta rápida. La Enfermedad X simboliza el potencial de una enfermedad inesperada y de rápida propagación, y resalta la necesidad de sistemas de salud flexibles y adaptables, así como capacidades de investigación para identificar, comprender y combatir patógenos desconocidos.
En esta práctica, va a aprender a estimar los rezagos epidemiológicos, el tiempo entre dos eventos epidemiológicos, utilizando un conjunto de datos simulado de la Enfermedad X.
La Enfermedad X es causada por un patógeno desconocido y se transmite directamente de persona a persona. Específicamente, la practica se centrará en estimar el período de incubación y el intervalo serial.
## 2. Agenda
Parte 1. Individual o en grupo.
Parte 2. En grupos de 6 personas. Construir estrategia de rastreo de contactos y aislamiento y preparar presentación de máximo 10 mins.
## 3. Conceptos claves
#### **3.1. Rezagos epidemiológicos: Período de incubación e intervalo serial**
En epidemiología, las distribuciones de rezagos se refieren a los *retrasos temporales* entre dos eventos clave durante un brote. Por ejemplo: el tiempo entre el inicio de los síntomas y el diagnóstico, el tiempo entre la aparición de síntomas y la muerte, entre muchos otros.
Este taller se enfocará en dos rezagos clave conocidos como el período de incubación y el intervalo serial. Ambos son cruciales para informar la respuesta de salud pública.
El [**período de incubación**]{.underline} es el tiempo entre la infección y la aparición de síntomas.
El [**intervalo serial**]{.underline} es el tiempo entre la aparición de síntomas entre el caso primario y secundario.
La relación entre estos parámetros tiene un impacto en si la enfermedad se transmite antes [(**transmisión pre-sintomática**)]{.underline} o después de que los síntomas [(**transmisión sintomática**)]{.underline} se hayan desarrollado en el caso primario (Figura 1).
Figura 1. Relación entre el período de incubación y el intervalo serial en el momento de la transmisión (*Adaptado de Nishiura et al. 2020)*
#### 3.2. Distribuciones comunes de rezagos y posibles sesgos
##### **3.2.1 Sesgos potenciales**
Cuando se estiman rezagos epidemiológicos, es importante considerar posibles sesgos:
[**Censura**]{.underline} significa que sabemos que un evento ocurrió, pero no sabemos exactamente cuándo sucedió. La mayoría de los datos epidemiológicos están "doblemente censurados" debido a la incertidumbre que rodea tanto los tiempos de eventos primarios como secundarios. No tener en cuenta la censura puede llevar a estimaciones sesgadas de la desviación estándar del resago (Park et al. en progreso).
[**Truncamiento a la derecha**]{.underline} es un tipo de sesgo de muestreo relacionado con el proceso de recolección de datos. Surge porque solo se pueden observar los casos que han sido reportados. No tener en cuenta el truncamiento a la derecha durante la fase de crecimiento de una epidemia puede llevar a una subestimación del rezago medio (Park et al. en progreso).
El sesgo [**dinámico (o de fase epidémica**)]{.underline} es otro tipo de sesgo de muestreo. Afecta a los datos retrospectivos y está relacionado con la fase de la epidemia: durante la fase de crecimiento exponencial, los casos que desarrollaron síntomas recientemente están sobrerrepresentados en los datos observados, mientras que durante la fase de declive, estos casos están subrepresentados, lo que lleva a la estimación de intervalos de retraso más cortos y más largos, respectivamente (Park et al. en progreso).
##### **3.2.2 Distribuciones de rezagos**
Tres distribuciones de probabilidad comunes utilizadas para caracterizar rezagos en epidemiología de enfermedades infecciosas (Tabla 1):
+----------------+-----------------------------------------+
| Distribución | Parámetros |
+================+:=======================================:+
| **Weibull** | `shape` y `scale` (forma y escala) |
+----------------+-----------------------------------------+
| **gamma** | `shape` y `scale` (forma y escala) |
+----------------+-----------------------------------------+
| **log normal** | `log mean` y `log standard deviation` |
| |(media y desviación estándar logarítmica)|
+----------------+-----------------------------------------+
: Tabla 1. Tres de las distribuciones de probabilidad más comunes para rezagos epidemiológicos.
## 4. Paquetes de *R* para la practica
En esta practica se usarán los siguientes paquetes de `R`:
- `dplyr` para manejo de datos
- `epicontacts` para visualizar los datos de rastreo de contactos
- `ggplot2` y `patchwork` para gráficar
- `incidence` para visualizar curvas epidemicas
- `rstan` para estimar el período de incubación
- `coarseDataTools` vía `EpiEstim` para estimar el intervalo serial
Instrucciones de instalación para los paquetes:
```{r install_packages, echo=FALSE}
## Packages on CRAN
#install.packages("dplyr")
#install.packages("incidence")
#install.packages("coarseDataTools")
#install.packages("EpiEstim")
#install.packages("ggplot2)
#install.packages("loo")
#install.packages("patchwork")
#install.packages("rstan")
#install.packages("parallel")
#
## Packages on GitHub
#remotes::install_github("reconhub/epicontacts")
```
Para cargar los paquetes, escriba:
```{r load_packages, warning=FALSE, message=FALSE}
library(dplyr)
library(epicontacts)
library(incidence)
library(coarseDataTools)
library(EpiEstim)
library(ggplot2)
library(loo)
library(patchwork)
library(rstan)
```
## 5. Datos
Esta práctica está partida en dos grupos para abordar dos enfermedades desconocidas con diferentes modos de transmisión.
Cargue los datos simulados que están guardados como un archivo .RDS, de acuerdo a su grupo asignado. Puede encontrar esta información en la carpeta [Enfermedad X](https://drive.google.com/drive/folders/1v8gMuEJx24ottM0VR4X_Ile60EV18ay3?usp=sharing). Descargue la carpeta, extráigala en el computador y abra el proyecto de R.
Hay dos elementos de interés:
- `linelist`, un archivo que contiene una lista de casos de la Enfermedad X, un caso por fila.
- `contacts`, un archivo con datos de rastreo de contactos que contiene información sobre pares de casos primarios y secundarios.
```{r load dat, echo=TRUE}
# Grupo 2
dat <- readRDS("data/practical_data_group2.RDS")
linelist <- dat$linelist
contacts <- dat$contacts
```
## 6. Exploración de los datos
#### **6.1. Exploración de los datos en `linelist`**
Empiece con `linelist`. Estos datos fueron recolectados como parte de la vigilancia epidemiológica rutinaria. Cada fila representa un caso de la Enfermedad X, y hay 7 variables:
- `id`: número único de identificación del caso
- `date_onset`: fecha de inicio de los síntomas del paciente
- `sex`: : M = masculino; F = femenino
- `age`: edad del paciente en años
- `exposure`: información sobre cómo el paciente podría haber estado expuesto
- `exposure_start`: primera fecha en que el paciente estuvo expuesto
- `exposure_end`: última fecha en que el paciente estuvo expuesto
::: {.alert .alert-secondary}
*💡 **Preguntas (1)***
- *¿Cuántos casos hay en los datos de `linelist`?*
- *¿Qué proporción de los casos son femeninos?*
- *¿Cuál es la distribución de edades de los casos?*
- *¿Qué tipo de información sobre la exposición está disponible?*
:::
```{r data exploration}
# Inspecionar los datos
head(linelist)
# P1
nrow(linelist)
# P2
table(linelist$sex)[2]/nrow(linelist)
# P3
summary(linelist$age)
# P4
table(linelist$exposure, exclude = F)[1]/nrow(linelist)
```
::: {.alert .alert-secondary}
*💡 **Discusión***
- ¿Por qué cree que falta la información de exposición en algunos casos?
- Ahora, grafique la curva epidémica. ¿En qué parte del brote cree que está (principio, medio, final)?
:::
```{r epi curve, warning = FALSE, message = FALSE}
i <- incidence(linelist$date_onset)
plot(i) +
theme_classic() +
scale_fill_manual(values = "purple") +
theme(legend.position = "none")
```
Parece que la epidemia todavía podría esta creciendo.
#### **6.2. Exploración de los datos de `rastreo de contactos`**
Ahora vea los datos de rastreo de contactos, que se obtuvieron a través de entrevistas a los pacientes. A los pacientes se les preguntó sobre actividades e interacciones recientes para identificar posibles fuentes de infección. Se emparejaron pares de casos primarios y secundarios si el caso secundario nombraba al caso primario como un contacto. Solo hay información de un subconjunto de los casos porque no todos los pacientes pudieron ser contactados para una entrevista.
Note que para este ejercicio, se asumirá que los casos secundarios solo tenían un posible infectante. En realidad, la posibilidad de que un caso tenga múltiples posibles infectantes necesita ser evaluada.
Los datos de rastreo de contactos tienen 4 variables:
- `primary_case_id`: número de identificación único para el caso primario (infectante)
- `secondary_case_id`: número de identificación único para el caso secundario (infectado)
- `primary_onset_date`: fecha de inicio de síntomas del caso primario
- `secondary_onset_date`: fecha de inicio de síntomas del caso secundario
```{r contacts}
x <- make_epicontacts(linelist = linelist,
contacts = contacts,
from = "primary_case_id",
to = "secondary_case_id",
directed = TRUE) # Esto indica que los contactos son directos (i.e., este gráfico traza una flecha desde los casos primarios a los secundarios)
plot(x)
```
::: {.alert .alert-secondary}
*💡 **Preguntas (2)***
- Describa los grupos (clusters).
- ¿Ve algún evento potencial de superpropagación (donde un caso propaga el patógeno a muchos otros casos)?
:::
# *`_______Pausa 1 _________`*
------------------------------------------------------------------------
## 7. Estimación del período de incubación
Ahora, enfoquese en el período de incubación. Se utilizará los datos del `linelist` para esta parte. Se necesitan ambos el tiempo de inicio de sintomas y el timpo de la posible exposición. Note que en los datos hay dos fechas de exposición, una de inicio y una de final. Algunas veces la fecha exacta de exposición es desconocida y en su lugar se obtiene la ventana de exposición durante la entrevista.
::: {.alert .alert-secondary}
*💡 **Preguntas (3)***
- ¿Para cuántos casos tiene datos tanto de la fecha de inicio de síntomas como de exposición?
- Calcule las ventanas de exposición. ¿Cuántos casos tienen una única fecha de exposición?
:::
```{r estimate ip, warning = FALSE}
ip <- filter(linelist, !is.na(exposure_start) &
!is.na(exposure_end))
nrow(ip)
ip$exposure_window <- as.numeric(ip$exposure_end - ip$exposure_start)
table(ip$exposure_window)
```
### 7.1. Estimación naive del período de incubación
Empiece calculando una estimación naive del período de incubación.
```{r naive}
# Máximo tiempo de período de incubación
ip$max_ip <- ip$date_onset - ip$exposure_start
summary(as.numeric(ip$max_ip))
# Mínimo tiempo de período de incubación
ip$min_ip <- ip$date_onset - ip$exposure_end
summary(as.numeric(ip$min_ip))
```
### 7.2. Censura estimada ajustada del período de incubación
Ahora, ajuste tres distribuciones de probabilidad a los datos del período de incubación teniendo en cuenta la censura doble. Adapte un código de `stan` que fue publicado por Miura et al. durante el brote global de mpox de 2022. Este método no tiene en cuenta el truncamiento a la derecha ni el sesgo dinámico.
Recuerde que el interés principal es considerar tres distribuciones de probabilidad: *Weibull*, *gamma* y *log normal* (Ver Tabla 1).
`Stan` es un programa de software que implementa el algoritmo Monte Carlo Hamiltoniano (HMC por su siglas en inglés de Hamiltonian Monte Carlo). HMC es un método de Monte Carlo de cadena de Markov (MCMC) para ajustar modelos complejos a datos utilizando estadísticas bayesianas.
#### **7.1.1. Corra el modelo en Stan**
Ajuste las tres distribuciones en este bloque de código.
```{r correct ip}
# Prepare los datos
earliest_exposure <- as.Date(min(ip$exposure_start))
ip <- ip |>
mutate(date_onset = as.numeric(date_onset - earliest_exposure),
exposure_start = as.numeric(exposure_start - earliest_exposure),
exposure_end = as.numeric(exposure_end - earliest_exposure)) |>
select(id, date_onset, exposure_start, exposure_end)
# Configure algunas opciones para ejecutar las cadenas MCMC en paralelo
# Ejecución de las cadenas MCMC en paralelo significa que se ejecutaran varias cadenas al mismo tiempo usando varios núcleos de su computador
options(mc.cores=parallel::detectCores())
input_data <- list(N = length(ip$exposure_start), # NNúmero de observaciones
tStartExposure = ip$exposure_start,
tEndExposure = ip$exposure_end,
tSymptomOnset = ip$date_onset)
# tres distribuciones de probabilidad
distributions <- c("weibull", "gamma", "lognormal")
# Código de Stan
code <- sprintf("
data{
int<lower=1> N;
vector[N] tStartExposure;
vector[N] tEndExposure;
vector[N] tSymptomOnset;
}
parameters{
real<lower=0> par[2];
vector<lower=0, upper=1>[N] uE; // Uniform value for sampling between start and end exposure
}
transformed parameters{
vector[N] tE; // infection moment
tE = tStartExposure + uE .* (tEndExposure - tStartExposure);
}
model{
// Contribution to likelihood of incubation period
target += %s_lpdf(tSymptomOnset - tE | par[1], par[2]);
}
generated quantities {
// likelihood for calculation of looIC
vector[N] log_lik;
for (i in 1:N) {
log_lik[i] = %s_lpdf(tSymptomOnset[i] - tE[i] | par[1], par[2]);
}
}
", distributions, distributions)
names(code) <- distributions
# La siguiente línea puede tomar ~1.5 min
models <- mapply(stan_model, model_code = code)
# Toma ~40 sec.
fit <- mapply(sampling, models, list(input_data),
iter=3000, # Número de iteraciones (largo de la cadena MCMC)
warmup=1000, # Número de muestras a descartar al inicio de MCMC
chain=4) # Número de cadenas MCMC a ejecutar
pos <- mapply(function(z) rstan::extract(z)$par, fit, SIMPLIFY=FALSE) # muestreo posterior
```
#### **7.1.2. Revisar si hay convergencia**
Ahora verifique la convergencia del modelo. Observe los valores de r-hat, los tamaños de muestra efectivos y las trazas MCMC. R-hat compara las estimaciones entre y dentro de cadenas para los parámetros del modelo; valores cercanos a 1 indican que las cadenas se han mezclado bien (Vehtari et al. 2021). El tamaño de muestra efectivo estima el número de muestras independientes después de tener en cuenta la dependencia en las cadenas MCMC (Lambert 2018). Para un modelo con 4 cadenas MCMC, se recomienda un tamaño total de muestra efectiva de al menos 400 (Vehtari et al. 2021).
Para cada modelo con distribución ajustada:
::: {.alert .alert-secondary}
*💡 **Preguntas (4)***
- ¿Los valores de r-hat son cercanos a 1?
- ¿Las 4 cadenas MCMC generalmente se superponen y permanecen alrededor de los mismos valores (se ven como orugas peludas)?
:::
#### **7.1.2.1. Convergencia para Gamma**
```{r convergence gamma}
print(fit$gamma, digits = 3, pars = c("par[1]","par[2]"))
rstan::traceplot(fit$gamma, pars = c("par[1]","par[2]"))
```
#### **7.1.2.2. Convergencia para log normal**
```{r convergence lognormal}
print(fit$lognormal, digits = 3, pars = c("par[1]","par[2]"))
rstan::traceplot(fit$lognormal, pars = c("par[1]","par[2]"))
```
#### **7.1.2.3. Convergencia para Weibull**
```{r convergence weibull}
print(fit$weibull, digits = 3, pars = c("par[1]","par[2]"))
rstan::traceplot(fit$weibull, pars = c("par[1]","par[2]"))
```
#### **7.1.3. Calcule los criterios de comparación de los modelos**
Calcule el criterio de información ampliamente aplicable (WAIC) y el criterio de información de dejar-uno-fuera (LOOIC) para comparar los ajustes de los modelos. El modelo con mejor ajuste es aquel con el WAIC o LOOIC más bajo. En esta sección también se resumirá las distribuciones y se hará algunos gráficos.
::: {.alert .alert-secondary}
*💡 **Preguntas (5)***
- ¿Qué modelo tiene mejor ajuste?
:::
####
```{r looic and waic ip, warning=FALSE, message=FALSE}
# Calcule WAIC para los tres modelos
waic <- mapply(function(z) waic(extract_log_lik(z))$estimates[3,], fit)
waic
# Para looic, se necesita proveer los tamaños de muestra relativos
# al llamar a loo. Este paso lleva a mejores estimados de los tamaños de
# muestra PSIS efectivos y del error de Monte Carlo
# Extraer la verosimilitud puntual logarítmica para la distribución Weibull
loglik <- extract_log_lik(fit$weibull, merge_chains = FALSE)
# Obtener los tamaños de muestra relativos efectivos
r_eff <- relative_eff(exp(loglik), cores = 2)
# Calcula LOOIC
loo_w <- loo(loglik, r_eff = r_eff, cores = 2)$estimates[3,]
# Imprimir los resultados
loo_w[1]
# Extraer la verosimilitud puntual logarítmica para la distribución gamma
loglik <- extract_log_lik(fit$gamma, merge_chains = FALSE)
r_eff <- relative_eff(exp(loglik), cores = 2)
loo_g <- loo(loglik, r_eff = r_eff, cores = 2)$estimates[3,]
loo_g[1]
# Extraer la verosimilitud puntual logarítmica para la distribución log normal
loglik <- extract_log_lik(fit$lognormal, merge_chains = FALSE)
r_eff <- relative_eff(exp(loglik), cores = 2)
loo_l <- loo(loglik, r_eff = r_eff, cores = 2)$estimates[3,]
loo_l[1]
```
#### **7.1.4. Reporte los resultados**
La cola derecha de la distribución del período de incubación es importante para diseñar estrategias de control (por ejemplo, cuarentena), los percentiles del 25 al 75 informan sobre el momento más probable en que podría ocurrir la aparición de síntomas, y la distribución completa puede usarse como una entrada en modelos matemáticos o estadísticos, como para pronósticos (Lessler et al. 2009).
Obtenga las estadísticas resumidas.
```{r reporting res}
# Necesitamos convertir los parámetros de las distribuciones a la media y desviación estándar del rezago
# En Stan, los parámetros de las distribuciones son:
# Weibull: forma y escala
# Gamma: forma e inversa de la escala (aka rate)
# Log Normal: mu y sigma
# Referencia: https://mc-stan.org/docs/2_21/functions-reference/positive-continuous-distributions.html
# Calcule las medias
means <- cbind(
pos$weibull[, 2] * gamma(1 + 1 / pos$weibull[, 1]), # media de Weibull
pos$gamma[, 1] / pos$gamma[, 2], # media de gamma
exp(pos$lognormal[, 1] + pos$lognormal[, 2]^2 / 2) # media de log normal
)
# Calcule las desviaciones estándar
standard_deviations <- cbind(
sqrt(pos$weibull[, 2]^2 * (gamma(1 + 2 / pos$weibull[, 1]) - (gamma(1 + 1 / pos$weibull[, 1]))^2)),
sqrt(pos$gamma[, 1] / (pos$gamma[, 2]^2)),
sqrt((exp(pos$lognormal[, 2]^2) - 1) * (exp(2 * pos$lognormal[, 1] + pos$lognormal[, 2]^2)))
)
# Imprimir los rezagos medios e intervalos creíbles del 95%
probs <- c(0.025, 0.5, 0.975)
res_means <- apply(means, 2, quantile, probs)
colnames(res_means) <- colnames(waic)
res_means
res_sds <- apply(standard_deviations, 2, quantile, probs)
colnames(res_sds) <- colnames(waic)
res_sds
# Informe la mediana e intervalos creíbles del 95% para los cuantiles de cada distribución
quantiles_to_report <- c(0.025, 0.05, 0.5, 0.95, 0.975, 0.99)
# Weibull
cens_w_percentiles <- sapply(quantiles_to_report, function(p) quantile(qweibull(p = p, shape = pos$weibull[,1], scale = pos$weibull[,2]), probs = probs))
colnames(cens_w_percentiles) <- quantiles_to_report
print(cens_w_percentiles)
# Gamma
cens_g_percentiles <- sapply(quantiles_to_report, function(p) quantile(qgamma(p = p, shape = pos$gamma[,1], rate = pos$gamma[,2]), probs = probs))
colnames(cens_g_percentiles) <- quantiles_to_report
print(cens_g_percentiles)
# Log normal
cens_ln_percentiles <- sapply(quantiles_to_report, function(p) quantile(qlnorm(p = p, meanlog = pos$lognormal[,1], sdlog= pos$lognormal[,2]), probs = probs))
colnames(cens_ln_percentiles) <- quantiles_to_report
print(cens_ln_percentiles)
```
Para cada modelo, encuentre estos elementos para el período de incubación estimado en la salida de arriba y escribalos abajo.
- Media e intervalo de credibilidad del 95%
- Desviación estándar e intervalo de credibilidad del 95%
- Percentiles (e.g., 2.5, 5, 25, 50, 75, 95, 97.5, 99)
- Los parámetros de las distribuciones ajustadas (e.g., shape y scale para distribución gamma)
#### **7.1.5. Grafique los resultados**
```{r plot ip}
# Prepare los resultados para graficarlos
df <- data.frame(
#Tome los valores de las medias para trazar la función de densidad acumulatica empirica
inc_day = ((input_data$tSymptomOnset-input_data$tEndExposure)+(input_data$tSymptomOnset-input_data$tStartExposure))/2
)
x_plot <- seq(0, 30, by=0.1) # Esto configura el rango del eje x (número de días)
Gam_plot <- as.data.frame(list(dose= x_plot,
pred= sapply(x_plot, function(q) quantile(pgamma(q = q, shape = pos$gamma[,1], rate = pos$gamma[,2]), probs = c(0.5))),
low = sapply(x_plot, function(q) quantile(pgamma(q = q, shape = pos$gamma[,1], rate = pos$gamma[,2]), probs = c(0.025))),
upp = sapply(x_plot, function(q) quantile(pgamma(q = q, shape = pos$gamma[,1], rate = pos$gamma[,2]), probs = c(0.975)))
))
Wei_plot <- as.data.frame(list(dose= x_plot,
pred= sapply(x_plot, function(q) quantile(pweibull(q = q, shape = pos$weibull[,1], scale = pos$weibull[,2]), probs = c(0.5))),
low = sapply(x_plot, function(q) quantile(pweibull(q = q, shape = pos$weibull[,1], scale = pos$weibull[,2]), probs = c(0.025))),
upp = sapply(x_plot, function(q) quantile(pweibull(q = q, shape = pos$weibull[,1], scale = pos$weibull[,2]), probs = c(0.975)))
))
ln_plot <- as.data.frame(list(dose= x_plot,
pred= sapply(x_plot, function(q) quantile(plnorm(q = q, meanlog = pos$lognormal[,1], sdlog= pos$lognormal[,2]), probs = c(0.5))),
low = sapply(x_plot, function(q) quantile(plnorm(q = q, meanlog = pos$lognormal[,1], sdlog= pos$lognormal[,2]), probs = c(0.025))),
upp = sapply(x_plot, function(q) quantile(plnorm(q = q, meanlog = pos$lognormal[,1], sdlog= pos$lognormal[,2]), probs = c(0.975)))
))
# Grafique las curvas de la distribución acumulada
gamma_ggplot <- ggplot(df, aes(x=inc_day)) +
stat_ecdf(geom = "step")+
xlim(c(0, 30))+
geom_line(data=Gam_plot, aes(x=x_plot, y=pred), color=RColorBrewer::brewer.pal(11, "RdBu")[11], linewidth=1) +
geom_ribbon(data=Gam_plot, aes(x=x_plot,ymin=low,ymax=upp), fill = RColorBrewer::brewer.pal(11, "RdBu")[11], alpha=0.1) +
theme_bw(base_size = 11)+
labs(x="Incubation period (days)", y = "Proportion")+
ggtitle("Gamma")
weibul_ggplot <- ggplot(df, aes(x=inc_day)) +
stat_ecdf(geom = "step")+
xlim(c(0, 30))+
geom_line(data=Wei_plot, aes(x=x_plot, y=pred), color=RColorBrewer::brewer.pal(11, "RdBu")[11], linewidth=1) +
geom_ribbon(data=Wei_plot, aes(x=x_plot,ymin=low,ymax=upp), fill = RColorBrewer::brewer.pal(11, "RdBu")[11], alpha=0.1) +
theme_bw(base_size = 11)+
labs(x="Incubation period (days)", y = "Proportion")+
ggtitle("Weibull")
lognorm_ggplot <- ggplot(df, aes(x=inc_day)) +
stat_ecdf(geom = "step")+
xlim(c(0, 30))+
geom_line(data=ln_plot, aes(x=x_plot, y=pred), color=RColorBrewer::brewer.pal(11, "RdBu")[11], linewidth=1) +
geom_ribbon(data=ln_plot, aes(x=x_plot,ymin=low,ymax=upp), fill = RColorBrewer::brewer.pal(11, "RdBu")[11], alpha=0.1) +
theme_bw(base_size = 11)+
labs(x="Incubation period (days)", y = "Proportion")+
ggtitle("Log normal")
(lognorm_ggplot|gamma_ggplot|weibul_ggplot) + plot_annotation(tag_levels = 'A')
```
En los gráficos anteriores, la línea negra es la distribución acumulativa empírica (los datos), mientras que la curva azul es la distribución de probabilidad ajustada con los intervalos de credibilidad del 95%. Asegúrese de que la curva azul esté sobre la línea negra.
::: {.alert .alert-secondary}
*💡 **Preguntas (6)***
- ¿Son los ajustes de las distribuciones lo que espera?
:::
####
# *`_______Pausa 2 _________`*
## 8. Estimación del intervalo serial
Ahora, estime el intervalo serial. Nuevamente, se realizará primero una estimación navie calculando la diferencia entre la fecha de inicio de síntomas entre el par de casos primario y secundario.
1. ¿Existen casos con intervalos seriales negativos en los datos (por ejemplo, el inicio de los síntomas en el caso secundario ocurrió antes del inicio de los síntomas en el caso primario)?
2. Informe la mediana del intervalo serial, así como el mínimo y el máximo.
3. Grafique la distribución del intervalo serial.
### 8.1. Estimación naive
```{r si naive}
contacts$diff <- as.numeric(contacts$secondary_onset_date - contacts$primary_onset_date)
summary(contacts$diff)
hist(contacts$diff, xlab = "Serial interval (days)", breaks = 25, main = "", col = "pink")
```
### 8.2. Estimación ajustada por censura
Ahora se estimará el intervalo serial utilizando una implementación del paquete `courseDataTools` dentro del paquete R `EpiEstim`. Este método tiene en cuenta la censura doble y permite comparar diferentes distribuciones de probabilidad, pero no se ajusta por truncamiento a la derecha o sesgo dinámico.
Se considerará tres distribuciones de probabilidad y deberá seleccionar la que mejor se ajuste a los datos utilizando WAIC o LOOIC. Recuerde que la distribución con mejor ajuste tendrá el WAIC o LOOIC más bajo.
Ten en cuenta que en `coarseDataTools`, los parámetros para las distribuciones son ligeramente diferentes que en rstan. Aquí, los parámetros para la distribución gamma son shape y scale (forma y escala) (<https://cran.r-project.org/web/packages/coarseDataTools/coarseDataTools.pdf>).
Solo se ejecutará una cadena MCMC para cada distribución en interés del tiempo, pero en la práctica debería ejecutar más de una cadena para asegurarse de que el MCMC converge en la distribución objetivo. Usará las distribuciones a priori predeterminadas, que se pueden encontrar en la documentación del paquete (ver 'detalles' para la función dic.fit.mcmc aquí:
(<https://cran.r-project.org/web/packages/coarseDataTools/coarseDataTools.pdf>).
#### 8.2.1. Preparación de los datos
```{r si epiestim}
# Formatee los datos de intervalos censurados del intervalo serial
# Cada línea representa un evento de transmisión
# EL/ER muestran el límite inferior/superior de la fecha de inicio de los síntomas en el caso primario (infector)
# SL/SR muestran lo mismo para el caso secundario (infectado)
# type tiene entradas 0 que corresponden a datos censurados doblemente por intervalo
# (ver Reich et al. Statist. Med. 2009)
si_data <- contacts |>
select(-primary_case_id, -secondary_case_id, -primary_onset_date, -secondary_onset_date,) |>
rename(SL = diff) |>
mutate(type = 0, EL = 0, ER = 1, SR = SL + 1) |>
select(EL, ER, SL, SR, type)
```
#### 8.2.2. Ajuste una distribución gamma para el SI
Primero, ajuste una distribución gamma al intervalo serial.
```{r gamma}
overall_seed <- 3 # semilla para el generador de números aleatorios para MCMC
MCMC_seed <- 007
# Ejecutaremos el modelo durante 4000 iteraciones con las primeras 1000 muestras descartadas como burning
n_mcmc_samples <- 3000 # número de muestras a extraer de la posterior (después del burning)
params = list(
dist = "G", # Ajuste de una distribución Gamma para el Intervalo Serial (SI)
method = "MCMC", # MCMC usando coarsedatatools
burnin = 1000, # número de muestras de burning (muestras descartadas al comienzo de MCMC)
n1 = 50, # n1 es el número de pares de media y desviación estándar del SI que se extraen
n2 = 50) # n2 es el tamaño de la muestra posterior extraída para cada par de media y desviación estándar del SI
mcmc_control <- make_mcmc_control(
seed = MCMC_seed,
burnin = params$burnin)
dist <- params$dist
config <- make_config(
list(
si_parametric_distr = dist,
mcmc_control = mcmc_control,
seed = overall_seed,
n1 = params$n1,
n2 = params$n2))
# Ajuste el SI
si_fit_gamma <- coarseDataTools::dic.fit.mcmc(
dat = si_data,
dist = dist,
init.pars = init_mcmc_params(si_data, dist),
burnin = mcmc_control$burnin,
n.samples = n_mcmc_samples,
seed = mcmc_control$seed)
```
Ahora observe los resultados.
```{r gamma check res}
# Verificar convergencia de las cadenas MCMC
converg_diag_gamma <- check_cdt_samples_convergence(si_fit_gamma@samples)
converg_diag_gamma
```
```{r gamma sample MCMC}
# Guardar las muestras MCMC en un dataframe
si_samples_gamma <- data.frame(
type = 'Symptom onset',
shape = si_fit_gamma@samples$var1,
scale = si_fit_gamma@samples$var2,
p50 = qgamma(
p = 0.5,
shape = si_fit_gamma@samples$var1,
scale = si_fit_gamma@samples$var2)) |>
mutate( # La ecuación de conversión se encuentra aquí: https://en.wikipedia.org/wiki/Gamma_distribution
mean = shape*scale,
sd = sqrt(shape*scale^2)
)
# Obtener la media, desviación estándar y 95% CrI
si_summary_gamma <-
si_samples_gamma %>%
summarise(
mean_mean = quantile(mean,probs=.5),
mean_l_ci = quantile(mean,probs=.025),
mean_u_ci = quantile(mean,probs=.975),
sd_mean = quantile(sd, probs=.5),
sd_l_ci = quantile(sd,probs=.025),
sd_u_ci = quantile(sd,probs=.975)
)
si_summary_gamma
```
```{r gamma check res si, echo=TRUE, results=FALSE}
# Obtenga las mismas estadísticas de resumen para los parámetros de la distribución
si_samples_gamma |>
summarise(
shape_mean = quantile(shape, probs=.5),
shape_l_ci = quantile(shape, probs=.025),
shape_u_ci = quantile(shape, probs=.975),
scale_mean = quantile(scale, probs=.5),
scale_l_ci = quantile(scale, probs=.025),
scale_u_ci = quantile(scale, probs=.975)
)
```
```{r gamma check res parms, results=FALSE}
# Necesita esto para hacer gráficos más tarde
gamma_shape <- si_fit_gamma@ests['shape',][1]
gamma_rate <- 1 / si_fit_gamma@ests['scale',][1]
```
#### 8.2.3. Ajuste de una distribución log normal para el intervalo serial
Ahora, ajuste una distribución log normal a los datos del intervalo serial.
```{r log normal}
# Ejecute el modelo durante 4000 iteraciones, descartando las primeras 1000 muestras como burning
n_mcmc_samples <- 3000 # número de muestras a extraer de la posterior (después del burning)
params = list(
dist = "L", # Ajustando una distribución log-normal para el Intervalo Serial (SI)
method = "MCMC", # MCMC usando coarsedatatools
burnin = 1000, # número de muestras de burning (muestras descartadas al comienzo de MCMC)
n1 = 50, # n1 es el número de pares de media y desviación estándar de SI que se extraen
n2 = 50) # n2 es el tamaño de la muestra posterior extraída para cada par de media y desviación estándar de SI
mcmc_control <- make_mcmc_control(
seed = MCMC_seed,
burnin = params$burnin)
dist <- params$dist
config <- make_config(
list(
si_parametric_distr = dist,
mcmc_control = mcmc_control,
seed = overall_seed,
n1 = params$n1,
n2 = params$n2))
# Ajuste del intervalo serial
si_fit_lnorm <- coarseDataTools::dic.fit.mcmc(
dat = si_data,
dist = dist,
init.pars = init_mcmc_params(si_data, dist),
burnin = mcmc_control$burnin,
n.samples = n_mcmc_samples,
seed = mcmc_control$seed)
```
Revise los resultados.
```{r log normal res conver}
# Revise la convergencia de las cadenas MCMC
converg_diag_lnorm <- check_cdt_samples_convergence(si_fit_lnorm@samples)
converg_diag_lnorm
```
```{r log normal res samples MCM, results=FALSE}
# Guarde las muestras de MCMC en un dataframe
si_samples_lnorm <- data.frame(
type = 'Symptom onset',
meanlog = si_fit_lnorm@samples$var1,
sdlog = si_fit_lnorm@samples$var2,