1. Introducción

Durante esta sesión vamos a explorar el comportamiento posiblemente dinámico de la Fuerza de Infección (FoI - Force of Infection) analizando encuestas serológicas poblacionales de prevalencia desagregadas por edad implementando modelos Bayesianos usando el paquete serofoi, una librería en lenguaje de programación R que permite implementar distintos modelos Bayesianos que estiman la FoI retrospectivamente.

2. Objetivos de aprendizaje

3. Conceptos básicos a desarrollar

En esta práctica se desarrollarán los siguientes conceptos:

3.1. Fuerza de Infección (FoI)

La Fuerza de Infección (FoI, por sus siglas en inglés), también conocida cómo la tasa de riesgo o la presión de infección, representa la tasa a la que los individuos susceptibles se infectan dado que estuvieron expuestos a un patógeno. En otras palabras, la FoI cuantifica el riesgo de que un individuo susceptible se infecte en un periodo de tiempo.

Como veremos a continuación, este concepto es clave en la modelamiento de enfermedades infecciosas. Para ilustrar el concepto, consideremos una población expuesta a un patógeno y llamemos \(P(t)\) a la proporción de individuos que han sido infectados al tiempo \(t\). Suponiendo que no hay sero-reversión, la cantidad \(1 - P(t)\) representa la cantidad de individuos susceptibles a la enfermedad, de tal forma que la velocidad de infección está dada por:

\[ \tag{1} \frac{dP(t)}{d t} = \lambda(t) (1- P(t)), \] en donde \(\lambda(t)\) representa la tasa a la que los individuos susceptibles se infectan, es decir la FoI, la cual se expresa como una tasa de infección por unidad de tiempo (días, meses, años, …). La ecuación diferencial (1) se asemeja a la de una reacción química en donde \(P(t)\) representa la proporción de sustrato que no ha entrado en contacto con un reactivo catalítico, por lo que este tipo de modelos se conocen como modelos catalíticos (Muench 1959).

Como las personas envejecen al mismo ritmo que pasa el tiempo, ecuaciones como la eq. (1) se pueden reescribir en términos de la edad \(a\) de los individuos de cada cohorte:

\[ \tag{2} \frac{dP(a)}{d a} = \lambda(a) (1- P(a)). \] Hacer este cambio nos habilita a usar, como base para estimar la FoI, encuestas serológicas que cumplan con los siguientes criterios de inclusión:

  • Basada en poblacionales (no hospitalaria).
  • Estudio de tipo transversal (un solo tiempo) .
  • Que indique la prueba diagnóstica utilizada.
  • Que identifique la edad del paciente en el momento de la encuesta.

En casos donde la FoI pueda considerarse constante, la ecuación (2) admite la siguiente solución analítica:

\[\tag{3} P(a) = 1 - \exp (-\lambda a), \]

en donde se tuvo en cuenta qué, al momento de la introducción del patógeno (\(t = 0\)), la proporción de individuos infectados fue \(P(0) = 0\) (condición inicial).

Figura 1. Curvas de prevalencia en función de la edad para distintos valores de FoI constante.

En este ejemplo hemos escogido, por simplicidad, que la FoI fuese constante; sin embargo, este no es necesariamente el caso. Identificar si la FoI sigue una tendencia constante o variante en el tiempo puede ser de vital importancia para la identificación y caracterización de la propagación de una enfermedad. Es acá donde el paquete de R serofoi juega un papel importante, ya que este permite estimar retrospectivamente la FoI de un patógeno, recuperando así su evolución temporal, por medio de modelos Bayesianos pre-establecidos.

Los modelos actuales del paquete serofoi asumen los siguientes supuestos biológicos:

  • No hay sero-reversión (sin pérdida de inmunidad).
  • La FoI no depende de la edad.
  • Bajos o nulos niveles de migración en las poblaciones.
  • Diferencias pequeñas entre las tasas de mortalidad de susceptibles e infectados.

3.2. Modelos Bayesianos

A diferencia del enfoque frecuentista, donde la probabilidad se asocia con la frecuencia relativa de la ocurrencia de los eventos, la estadística Bayesiana se basa en el uso de la probabilidad condicional de los eventos respecto al conocimiento (o estado de información) que podamos tener sobre los datos o sobre los parámetros que queramos estimar.

Por lo general, cuando proponemos un modelo lo que buscamos es disminuir la incertidumbre sobre algún parámetro, de tal forma que nos aproximemos a su valor tan óptimamente como nos lo permita nuestro conocimiento previo del fenómeno y de las mediciones realizadas (datos).

La inferencia Bayesiana se sustenta en el teorema de Bayes, el cual establece que: dado un conjunto de datos \(\vec{y} = (y_1, …, y_N)\), el cual representa un único evento, y la variable aleatoria \(\theta\), que representa un parámetro de interés para nosotros, la distribución de probabilidad conjunta de las variables aleatorias asociadas está dada por:

\[ \tag{4} p(\vec{y}, \theta) = p(\vec{y} | \theta) p(\theta) = p(\theta | \vec{y}) p(\vec{y}), \] de donde se desprende la distribución aposteriori de \(\theta\), es decir una versión actualizada de la distribución de probabilidad de la FoI condicionada a nuestros datos:

\[\tag{5} p(\theta, \vec{y}) = \frac{p(\vec{y} | \theta) p(\theta)}{p(\vec{y})}, \] La distribución \(p(\vec{y} | \theta)\), que corresponde a la información interna a los datos condicionada al valor del parámetro \(\theta\), suele estar determinada por la naturaleza del experimento: no es lo mismo escoger pelotas dentro de una caja reemplazándolas que dejándolas por fuera de la caja. En el caso particular de la FoI, contamos con datos como el número total de encuestas por edad y el número de casos positivos, por lo que es razonable asignar una distribución binomial a la probabilidad, como veremos a continuación.

3.3.1. Modelo FoI constante

En nuestro caso particular, el parámetro que queremos estimar es la FoI (\(\lambda\)). La distribución de probabilidad apriori de \(\lambda\) representa nuestras suposiciones informadas o el conocimiento previo que tengamos sobre el comportamiento de la FoI. En este contexto, el estado de mínima información sobre \(\lambda\) está representado por una distribución uniforme: \[\tag{6} \lambda \sim U(0, 2), \] lo que significa que partimos de la premisa de que todos los valores de la fuerza de infección entre \(0\) y \(2\) son igualmente probables. Por otro lado, de la teoría de modelos sero-catalíticos sabemos que la seroprevalencia en un año dado está descrita por un proceso cumulativo con la edad (Hens et al. 2010):

\[\tag{7} P(a, t) = 1 - \exp\left( -\sum_{i=t-a+1}^{t} \lambda_i \right), \] en donde \(\lambda_i\) corresponde a la FoI al tiempo \(t\). Como en este caso la FoI es constante, la ec. (7) se reduce a:

\[\tag{8} P(a, t) = 1 - \exp\left( -\lambda a \right), \] Si \(n(a, t_{sur})\) es el número de muestras positivas por edad obtenidas en un estudio serológico desarrollado en el año \(t_{sur}\), entonces podemos estimar la distribución de los casos seropositivos por edad condicionada al valor de \(\lambda\) como una distribución binomial: \[\tag{9} \lambda(t)\sim normal(\lambda(t-1), \sigma) \\ \lambda(t=1) \sim normal(0, 1) \]

3.3.2. Modelo FOI dependientes del tiempo

Actualmente, serofoi permite la implementación de dos modelos dependientes del tiempo: uno de variación lenta de la FoI (tv-normal) y otro de variación rápida (tv-normal-log) de la FoI.

Cada uno de ellos se basa en distintas distribuciones previas para \(\lambda\), las cuales se muestran en la tabla 1.

Model Option Probability of positive case at age \(a\) Prior distribution
constant \(\sim binom(n(a,t), P(a,t))\) \(\lambda\sim uniform(0,2)\)
tv_normal \(\sim binom(n(a,t), P(a,t))\) \(\lambda\sim normal(\lambda(t-1),\sigma)\\ \lambda(t=1)\sim normal(0,1)\)
tv_normal_log \(\sim binom(n(a,t), P(a,t))\) \(\lambda\sim normal(log(\lambda(t-1)),\sigma)\\ \lambda(t=1)\sim normal(-6,4)\)

Tabla 1. Distribuciones a priori de los distintos modelos soportados por serofoi. \(\sigma\) representa la desviación estándar.

Como se puede observar, las distribuciones previas de \(\lambda\) en ambos casos están dadas por distribuciones Gaussianas con desviación estándar \(\sigma\) y centradas en \(\lambda\) (modelo de variación lenta - tv_normal) y \(\log(\lambda)\) (modelo de variación rápida - tv_normal_log). De esta manera, la FoI en un tiempo \(t\) está distribuida de acuerdo a una distribución normal alrededor del valor que esta tenía en el tiempo inmediatamente anterior. El logaritmo en el modelo tv_normal_log permite identificar cambios drásticos en la tendencia temporal de la FoI.

4. Instalación de serofoi

Para llevar a cabo los análisis del taller es necesario instalar y cargar el paquete serofoi. Previo a la instalación de serofoi, cree un proyecto en R en la carpeta de su escogencia en su máquina local; esto con el fin de organizar el espacio de trabajo en donde se guardarán los códigos que desarrolle durante la sesión. Para ello, siga los pasos mostrados en el material de apoyo Introducción a R y RStudio. Una vez creado el proyecto, abra el archivo de extensión .Rproj en RStudio y ejecute las siguientes líneas de código:

if(!require("remotes")) install.packages("remotes")
remotes::install_github("epiverse-trace/serofoi")

5. Caso de uso 1: Enfermedad de Chagas

La enfermedad de Chagas es una infección parasitaria endémica de Latinoamérica que es causada por el protozoario Trypanosoma cruzi. Ésta es transmitida al humano por medio de las heces de insectos triatominos infectados, los cuales han estado presentes en Latinoamérica por miles de años. Los gusanos triatominos han desarrollado hábitos domiciliarios, cohabitando viviendas con humanos. El uso de insecticida es una de las estrategias de control primarias contra la enfermedad de Chagas, dado que reduce efectivamente la población de insectos triatominos en ambientes domésticos. Durante esta sesión estudiaremos una encuesta serológica recolectada en 2012, la cual midió la seroprevalencia de anticuerpos de tipo IgG contra la infección de Trypanosoma cruzi en una zona rural de Colombia. (Cucunubá et al. 2017) Una vez cargada la librería, ejecute el siguiente código para cargar y ver los datos de la encuesta:

data(chagas2012)
str(chagas2012)
## 'data.frame':    72 obs. of  9 variables:
##  $ survey  : chr  "COL-035-93" "COL-035-93" "COL-035-93" "COL-035-93" ...
##  $ total   : num  34 25 35 29 36 23 34 18 21 24 ...
##  $ counts  : num  0 0 1 0 0 0 0 0 0 0 ...
##  $ age_min : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ age_max : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ tsur    : num  2012 2012 2012 2012 2012 ...
##  $ country : chr  "COL" "COL" "COL" "COL" ...
##  $ test    : chr  "ELISA" "ELISA" "ELISA" "ELISA" ...
##  $ antibody: chr  "IgG anti-T.cruzi" "IgG anti-T.cruzi" "IgG anti-T.cruzi" "IgG anti-T.cruzi" ...

Para visualizar gráficamente los datos podemos utilizar la función plot_seroprev():

plot_seroprev(serodata = chagas2012)

Este comando nos muestra la gráfica de seropositividad de la encuesta. El tamaño de los puntos indica el número total de pruebas para cada edad, mientras que la altura representa la proporción de dichas pruebas que son positivas junto con su respectiva barra de error binomial.

¿Qué puede decir de estos datos? ¿Qué comportamiento espera de la FoI en este caso?


5.2. Estimación retrospectiva de la FoI

Como la enfermedad de Chagas es endémica, solo usaremos los modelos endémicos de serofoi en este caso, es decir constant y tv_normal. Antes de implementar los modelos, es necesario preparar nuestros datos para el proceso de modelamiento. Para ello, usamos la función prepare_serodata, la cual añade a la encuesta el marcador de edad de cada cohorte (age_mean_f), su respectivo año de nacimiento (birth_year), el tamaño de muestra (sample_size) y la prevalencia observada junto con su respectivo intervalo de confianza (prev_obs, prev_obs_lower, prev_obs_upper):

chagas2012p <- prepare_serodata(chagas2012)

Los modelos se pueden implementar haciendo uso de la función run_seromodel(), la cual permite especificar el nombre del modelo a utilizar (foi_model) y el número de iteraciones del algoritmo que implementa el modelo; entre más alto sea el número de iteraciones, mejor será la estimación de la FoI (a costas de incrementar el tiempo de implementación). En este caso, usaremos 800 iteraciones e implementaremos el modelo constante:

chagas_constant <- run_seromodel(serodata = chagas2012p, 
                                 foi_model = 'constant',
                                 n_iters = 800)

Es posible obtener un sumario de los resultados del modelo por medio de la función extract_seromodel_summary():

chagas_constant_summary <- extract_seromodel_summary(chagas_constant, 
    serodata = chagas2012p)
print(t(chagas_constant_summary))
##           [,1]              
## foi_model "constant"        
## dataset   "COL-035-93"      
## country   "COL"             
## year      "2012"            
## test      "ELISA"           
## antibody  "IgG anti-T.cruzi"
## n_sample  "747"             
## n_agec    "72"              
## n_iter    "800"             
## elpd      "-92.66"          
## se        "6.41"            
## converged "Yes"

5.3. Visualización de los resultados

Para visualizar los resultados de la implementación de una forma compacta podemos usar la función plot_seromodel, la cual grafica el ajuste de la seroprevalencia, la estimación retrospectiva de la FoI y un criterio de convergencia del modelo:

chagas_constant_plot <- plot_seromodel(chagas_constant, 
                                       serodata = chagas2012p,
                                       size_text = 6)
print(chagas_constant_plot)

Ejercicio:

Repita estos mismos pasos para el modelo de variación lenta de la FoI y responda las siguientes preguntas:

  • ¿En qué se diferencian los ajustes de la seroprevalencia entre los dos casos?




  • Gráficamente, ¿Cuál de los dos modelos parece ajustarse mejor a estos datos?


Pista: Debería obtener la siguiente gráfica para el modelo de variación lenta:

Dos criterios estadísticos que se pueden usar para determinar qué modelo se ajusta mejor a los datos son los siguientes:

  • Menor desviación estándar: A menor desviación estándar, menor es la dispersión de los datos, lo que sugiere un grado más alto de precisión en la estimación.
  • Mayor valor de elpd: El elpd (expected log-pointwise predictive density) es una medida de la capacidad predictiva del modelo; en cuanto mayor sea su valor (menos negativo), mejor será su precisión predictiva.

Según estos criterios, y teniendo en cuenta el contexto de la encuesta:

  • ¿Cuál modelo se ajusta mejor a los datos de Chagas?




  • ¿Tuvo algún efecto la intervención realizada con insecticidas en la zona?




6. Caso de uso 2: Chikungunya

Chikungunya es una enfermedad viral descubierta en Tanzania en el año 1952 que se caracteriza por causar fuertes fiebres, dolor articular, jaqueca, sarpullidos, entre otros síntomas. Desde el año 2013, casos de esta enfermedad comenzaron a ser reportados en América y desde entonces la enfermedad ha sido endémica en varios países latinoamericanos.

En esta sección vamos a analizar una encuesta serológica realizada entre octubre y diciembre de 2015 en Bahía, Brasil, la cual fue realizada poco tiempo después de una epidemia de esta enfermedad en la zona. Nuestro objetivo es caracterizar la propagación de la enfermedad por medio de la implementación de los distintos modelos y determinar cuál de estos describe mejor la situación. Primero, carguemos y preparemos los datos que utilizaremos en este análisis. La base chik2015 contiene los datos correspondientes a esta encuesta serológica:

data(chik2015)
chik2015p <- prepare_serodata(chik2015)
str(chik2015p)
## 'data.frame':    4 obs. of  16 variables:
##  $ survey        : chr  "BRA 2015(S019)" "BRA 2015(S019)" "BRA 2015(S019)" "BRA 2015(S019)"
##  $ total         : int  45 109 144 148
##  $ counts        : num  17 55 63 69
##  $ age_min       : num  1 20 40 60
##  $ age_max       : num  19 39 59 79
##  $ tsur          : num  2015 2015 2015 2015
##  $ country       : Factor w/ 256 levels "ABW","AFG","AGO",..: 33 33 33 33
##  $ test          : chr  "ELISA" "ELISA" "ELISA" "ELISA"
##  $ antibody      : chr  "IgG anti-CHIKV" "IgG anti-CHIKV" "IgG anti-CHIKV" "IgG anti-CHIKV"
##  $ reference     : chr  "Dias et al. 2018" "Dias et al. 2018" "Dias et al. 2018" "Dias et al. 2018"
##  $ age_mean_f    : num  10 29 49 69
##  $ sample_size   : int  446 446 446 446
##  $ birth_year    : num  2005 1986 1966 1946
##  $ prev_obs      : num  0.378 0.505 0.438 0.466
##  $ prev_obs_lower: num  0.238 0.407 0.355 0.384
##  $ prev_obs_upper: num  0.535 0.602 0.523 0.55

Para correr el modelo de FoI constante y visualizar los resultados del mismo, corra las siguientes líneas de código:

chik_constant <- run_seromodel(serodata = chik2015p,
                               foi_model = "constant",
                               n_iters = 1000)

chik_constant_plot <- plot_seromodel(chik_constant, 
                                     serodata = chik2015p, 
                                     size_text = 12)

Ahora, corra los modelos dependientes del tiempo con n_iters =1500. Luego visualice conjuntamente las tres gráficas por medio de la función plot_grid() del paquete cowplot:

install.packages("cowplot")

cowplot::plot_grid(chik_constant_plot, 
                   chik_normal_plot, 
                   chik_normal_log_plot,
                   ncol = 3)

NOTA: Debido a que el número de trayectorias es relativamente alto, con el fin de asegurar la convergencia de los modelos, el tiempo de cómputo de los modelos dependientes del tiempo puede tardar varios minutos.

Pista: Debería obtener la siguiente gráfica:

Según los criterios usados anteriormente, responda:

¿Cuál de los tres modelos se ajusta mejor a esta encuesta serológica?




¿Cómo interpreta estos resultados?


Contribuciones

  • Nicolás T. Domínguez
  • Zulma M. Cucunuba
  • Geraldine Sarai Gómez: Sugerencias

Contribuciones son bienvenidas vía pull requests.

Referencias


Cucunubá, Zulma M, Pierre Nouvellet, Lesong Conteh, Mauricio Javier Vera, Victor Manuel Angulo, Juan Carlos Dib, Gabriel Jaime Parra -Henao, and Marı́a Gloria Basáñez. 2017. “Modelling Historical Changes in the Force-of-Infection of Chagas Disease to Inform Control and Elimination Programmes: Application in Colombia.” BMJ Global Health 2 (3). https://doi.org/10.1136/bmjgh-2017-000345.
Hens, Niel, Marc Aerts, Christel Faes, Ziv Shkedy, O Lejeune, P Van Damme, and P Beutels. 2010. “Seventy-Five Years of Estimating the Force of Infection from Current Status Data.” Epidemiology & Infection 138 (6): 802–12.
Muench, Hugo. 1959. Catalytic Models in Epidemiology. Harvard University Press.