Introducción

El Dengue es la infección viral transmitida por artrópodos más importante de las últimas décadas para los seres humanos, tanto así que se estima que cerca de 2500 millones de personas están en riesgo de desarrollar dicha enfermedad (^{Chien, Lung-Chang. y cols., 2014}), el riesgo aumenta en las regiones tropicales o subtropicales, ya que el mosquito Aedes aegypti, su principal vector de transmisión, es una especie endémica para la mayoría de estas regiones. A pesar de que existen modelos teóricos cuya finalidad es ajustar un modelo de propagación, los cuales dan gran importancia a variables como la temperatura y la precipitación, y buscan determinar patrones de transferencia; todavía hace falta mayor evidencia empírica para lograr mejores conclusiones.

Para el desarrollo de sistemas de alerta tempranos, que ayuden a enfrentar la transmisión de la enfermedad, es esencial entender las relaciones empíricas entre algunos factores meteorológicos, geográficos, socioeconómicos y la fiebre del dengue (^{Chien, Lung-Chang. Y cols., 2014}). En este artículo se usa un modelo espacial Condicional Autoregresivo (CAR, por sus siglas en inglés) para estudiar el riesgo relativo en cada municipio de la Isla Principal de Puerto Rico. En el modelo propuesto analiza la relación de variables como: temperatura, precipitación y altitud, con el riesgo relativo de transmisión del dengue, aplicado a datos recolectados en el año 2014.

Datos utilizados en la investigación

El modelo construido incorpora el riesgo relativo de contraer la demás e del dengue por municipio de Puerto Rico, este riesgo se define como:

representa el riesgo relativo para las personas que habitan en el municipio i-ésimo, por ejemplo los datos de la semana 1 de 2014 arrojan el comportamiento descrito en la Figura 1 que se presenta a continuación:

Figura 1: Mapa de colores para el riesgo relativo en semana 1 

es el número de presuntos casos de dengue para el municipio i. Este dato se obtuvo de los reportes semanales de casos presuntos de Dengue, reportes elaborados en conjunto por el Centro para el Control y Prevención de Enfermedades de Estados Unidos (CDC, por sus siglas en ingles) y el Departamento de Salud de Puerto Rico. Específicamente, fueron obtenidos en la página de internet de dicho departamento de salud, ver (CDC, Subdivisión de Dengue y Departamento de Salud, PR, 2014),

es la población del municipio i, según el censo del año 2010, disponibles en la página de la Junta de Planificación del Gobierno de Puerto Rico (^{Bureau, 2011}),

es el número de casos presuntos totales de Puerto Rico, es decir que es la suma del número presunto de casos de todos los municipios,

es la población total en Puerto Rico según el censo del año 2010.

Las variables independientes (también llamadas covariables) utilizadas dentro del modelo son:

Modelación con datos de área

En el presente apartado se estudian 5 modelos distintos, en el sentido que ajustan los datos con enfoques basados en diferentes supuestos, cada modelo realiza un ajuste para los casos de dengue usando regresión de Poisson (ya que se modelan datos de conteo), sin embargo, dos no son espaciales, y los resultados que se dieron en algunas semanas de estudio nos indican que no existe necesidad de modelos espaciales para estos datos. Los modelos son:

Datos y análisis exploratorio: La región de estudio es la isla principal de Puerto Rico, que consta de 76 municipios, los cuales tienen un promedio de 48877 habitantes según el censo de 2010. Las variables que se usan en cada modelo son: altitud, nivel o porcentaje de pobreza, número de habitantes (población), superficie por municipio y densidad poblacional, dichas variables no varían entre semanas. Los datos se resumen en el Tabla 2, mostrando los percentiles de sus distribuciones.

Por otro lado están: precipitación medida en decimas de milímetros, temperatura mínima medida en grados Celsius, el índice de vegetación mejorado (EVI) y los casos presuntos de dengue, que se usan en lugar de los casos confirmados, ya que estos están fuertemente correlacionados y además son mas sensitivos porque aproximadamente el 60% de los casos confirmados tienen muestras poco adecuadas para un diagnostico definitivo (^{Johansson, Michael A. y cols., 2009}) y el riesgo relativo, este último calculado como el cociente de los casos presuntos y el riesgo de la población total de cada municipio.

La variable de respuesta en el presente estudio es dada por los casos de dengue, la cual se asume que sigue una distribución de Poisson, es decir, , tal que , en la última ecuación se cumple que y este término es un offset (número esperado de casos en el municipio k) usado para controlar el tamaño de la población, mientras que el termino contiene el efecto de las covariables, un efecto aleatorio espacialmente estructurado (en el caso de los modelos espaciales) y un efecto aleatorio sin estructura espacial.

En la Tabla 3 se presenta un resumen para tres de las semanas de estudio en el presente trabajo, donde se muestra que en los datos de casos hay cierta estabilidad entre los municipios, ya que al menos un 75 % de los mismos presentan un máximo de 2 casos de dengue en las tres semanas, y una media entre 1 y 2 casos. En cuanto a las variables de precipitación, temperatura e índice de vegetación debemos recordar que son el resultado de interpolación espacial mientras que el índice de vegetación se obtuvo por interpolación de la serie temporal, ver (^{Hernández-González, 2017}), por lo que para los modelos ajustados, estas variables podrían aportar un error aún mayor con respecto a si fueran variables medidas directamente (muestras).

En la Figura 1 se puede ver existe un municipio (Patillas) con un riesgo muy superior, también sus vecinos, esto se debe a que dichos municipios según los informes del Departamento de Salud, pertenecieron a un proyecto de vigilancia aumentada (^{CDC, Sub- división de Dengue y Departamento de Salud, PR, 2014}). Para las primeras 3 semanas el municipio de Patillas presenta el mayor número de casos de dengue, San Juan es en esas semanas posee el segundo lugar con mayor número de casos, pero la población de San Juan es poco mas de 20 veces la población de Patillas lo que se traduce en un riesgo relativo muy alto para el este municipio. Dado que el objetivo principal es ajustar modelos que usen todos los datos de la isla principal de Puerto Rico, y teniendo en cuenta la información del Departamento de Salud de Puerto Rico, la cual sugiere la existencia de mediciones muy diferentes para algunos municipios, se decidió usar la información completa de la isla, sin hacer ninguna valoración en el uso de herramientas para detectar valores atípicos (“outliers”) en los datos del estudio. Este comportamiento también se presento en las demás semanas de estudio.

Un tema importante a tratar es el de la autocorrelación espacial de las variables en estudio, en particular para el riesgo relativo, ya que es la variable dependiente del modelo, para llevar a cabo dicha exploración se usan pruebas de hipótesis para los índices de Moran y Geary. Los supuestos subyacentes en las pruebas son sensitivos a la matriz de peso utilizada, por ejemplo, una matriz de pesos estandarizada por las incrementa la influencia en los vínculos para observaciones con pocos vecinos, mientras que una matriz binaria varía la influencia de las observaciones: las que tienen muchos vecinos están sobre-ponderadas en comparación con aquellas con pocos vecinos (^{Bivand, Roger S. y cols., 2013}); las pruebas tienen como hipótesis nula que no existe autocorrelación espacial y de hipótesis alternativa que existe autocorrelación espacial, ya sea positiva o negativa, ambas pruebas se hicieron usando las funciones moran.mc() y geary.mc() del paquete spdep (Bivand, Roger. y cols., 2015, ²⁰¹³) de R-CRAN. Los resultados de dichas pruebas se presentan en la Tabla 4.

La prueba para altitud da como resultado un valor de 0,3213 para el estadístico con un p-valor de 0,0003, lo que indica que hay evidencia para rechazar la hipótesis nula de no autocorrelación espacial en dicha variable, por otro lado para el registro de pobreza se tiene un índice de Moran de 0,5684, teniendo en este caso un p-valor de 0,0001 por lo que se rechaza la hipótesis nula también es este caso. En la Tabla 4 se puede ver que con la excepción de las semanas 2 y 4 para el riesgo relativo, semana 5 para precipitación, las pruebas indican que las variables están autocorrelacionadas espacialmente. Es importante aclarar que dichas pruebas se realizaron con una matriz de pesos estandarizada por filas, sin embargo, también se hizo el ejercicio con la matriz binaria y los resultados obtenidos son similares.

Modelos no espaciales: Los modelos ajustados sin tomar en cuenta una estructura espacial son dos: modelo lineal generalizado (GLM), el cual se ajusto usando máxima verosimilitud y por medio de cadenas de Markov de Monte Carlo, y un modelo lineal generalizado mixto (GLMM), también llamado Independiente en este trabajo (^{Lee, 2013}), ajustado con MCMC. La forma que toma cada uno viene dada por:

Donde en el GLM se tiene que:

;

Mientras que el modelo Independiente se rige por la ecuación:

donde es el offset que controla el tamaño de la población, es el intercepto, es el coeficiente para la variable de altitud, es el coeficiente para el porcentaje de pobreza, es la proporción del índice de vegetación en el modelo, es un coeficiente relacionado a la densidad poblacional, es el coeficiente para la precipitación promedio en el municipio, es el coeficiente para la temperatura mínima y es una variable aleatoria con distribucion normal de media cero y varianza , a dicha varianza se le asigna una previa U(0,1000), lo anterior con el objetivo de que sea una previa poco informativa en el modelo.

En el primer escenario se revisa el ajuste realizado al modelo (2) usando máxima verosimilitud y MCMC, donde es importante destacar que la diferencia entre los estimadores es muy poca, sin embargo, con máxima verosimilitud no se puede hablar de un intervalo de predicción para cada uno de los estimadores. Para realizar el ajuste con máxima verosimilitud se uso la función glm() del paquete “stats" (R Core Team, 2016), mientras que el ajuste bayesiano se realizo con la función stan_glm( ), la cual pertenece al paquete “rstanarm" (^{Gabry y Goodrich, 2016}); las previas utilizadas por dicha función para los parámetros de ajuste asociados a las covariables son distribuciones normales N(0;100).

La primera conclusión que se puede dar, es que las covariables de EVI y Densidad no aportan mayor significancia al modelo, ya que para las 8 semanas de estudio el valor de la mediana (y también la media) de los parámetros correspondientes eran cero o muy cercanos a cero, además que el intervalo de predicción para dichos valores es muy pequeño y siempre contiene al 0. Los parámetros para las covariables EVI y densidad se mantienen muy cerca de cero o son cero para las demás semanas de estudio, aún así se hace el ajuste de los otros modelos teniendo en cuenta esas variables, pero el resultado obtenido en cada uno los ajustes fue el mismo, por lo que se volvió a ajustar los datos de todas las semanas con todos los modelos propuestos en la Tabla 1, pero sin dichas covariables, es decir, los nuevos modelos a ajustar son:

;

En tanto que el modelo Independiente es:

;

Donde cada elemento en las ecuaciones anteriores ya se han explicado. En el presente estudio los intervalos de predicción al 95 % del modelo lineal generalizado, sugieren que las covariables de mayor peso son altitud y pobreza ya que para todas las semanas ninguno contiene al cero, mientras que para precipitación y temperatura mínima el cero estaba presente en cada intervalo de predicción. Por lo anterior, en caso de existir una estructura espacial que no está siendo tomada en cuenta en dicho modelo, altitud y pobreza tienen un rol importante en explicar el modelo espacial en el riesgo relativo de contraer dengue en alguno de los municipios de Puerto Rico(^{Lee, 2013}). En la Tabla 5 se presentan tales intervalos.

Con el nuevo ajuste al parecer para estos datos el GLM no es una buena elección según el criterio de información de la devianza, ya que es menor en el modelo independiente para todas las semanas, por otro lado, el WAIC sugiere que las semanas 2 y 3 se ajustan mejor con un GLM en lugar de un GLMM, y las restantes semanas con el Independiente (en adelante GLMM).

A partir de aquí, se hace la comparación y selección de modelos usando únicamente el índice de Watanabe-Akaike, porque el DIC (^{Spiegelhalter, David J., y cols, 2002}) no tiene un enfoque bayesiano en su totalidad ya que se basa en un estimador puntual, mientras que el WAIC tiene un enfoque completamente bayesiano en el sentido que usa la distribución posterior completa, siendo además una mejora para el DIC en modelos bayesianos y es asintóticamente equivalente a la validación cruzada propuesta por ^{Vehtari en (Vehtari, Aki. y cols., 2016}), conocida como loo-cv (Leave-One-Out Cross Validation). Otras de las desventajas del DIC es que puede producir estimadores negativos del número efectivo de parámetros en un modelo y no está definido para modelos singulares. Finalmente, a diferencia del DIC, el WAIC es invariante a la parametrización y también es útil en la comparación para modelos con distribuciones singulares (Vehtari, Aki. y cols., 2016).

En el GLM, los estimadores se mantienen estables en las semanas de estudio, se puede ver en la Tabla 6 que el valor que toma para altitud es prácticamente constante (con 4 dígitos decimales), y los valores asociados a las demás covariables no tienen cambios grandes de semana a semana. Otro resultado obtenido desde el modelo lineal generalizado es que riesgo relativo es inversamente relacionado con altitud, ya que los parámetros de ajuste son negativos para las 8 semanas, eso quiere decir que a mayor altitud existe menor riesgo de contraer dengue, resultado que también obtienen otros autores, por ejemplo (^{Mena, Nelson. y cols., 2011}). Según los resultados del GLM, el índice de pobreza está relacionado de manera inversa al riesgo de contraer el virus, sin embargo, en el modelo Independiente se concluye que a pesar de la fuerte asociación que hay entre riesgo relativo y pobreza, no existe un efecto consistente entre las semanas, porque en algunas semanas el efecto es negativo y en otras es positivo, dicho resultado es semejante al que obtuvieron en un estudio del 2009 realizado por ^{Johansson, Dominici y Glass en (Johansson, Michael A. y cols., 2009}).

Resultados de modelos espaciales: Se ajustan 3 modelos espaciales: intrínseco, Besag y Leroux. El primero asume que existe completa autocorrelación espacial en la variable dependiente, el siguiente es una suma entre el modelo intrínseco y el independiente (5) y el último se puede ver como una ponderación entre los modelos intrínseco e independiente. La variable dependiente es el número de casos de dengue para cada municipio, la que suponemos que sigue una distribución de Poisson, es decir:

, donde

Para los modelos intrínseco y Leroux se cumple ( ^{Lee, 2013} ):

Intrínseco:

Leroux:

Los tres modelos se ajustan con funciones del paquete CARBayes (^{Lee, 2013}), presente en el programa estadístico R-CRAN (^{R Core Team, 2016}). Se realizan 50000 simulaciones con Monte Carlo vía cadenas de Markov, de las que se queman 10000 con el fin de haber alcanzado estabilidad o equilibrio en la cadena, donde los valores de los parámetros de la distribución gamma inversa son los preestablecidos en el paquete, y la matriz W utilizada es la matriz binaria definida como wij = 1 si los municipios ij comparten alguna fracción de sus fronteras y wij = 0 en otro caso. Además, para explorar convergencia en las cadenas de Markov se usa el diagnóstico de Geweke, con el que se concluye que en todas las semanas cada parámetro ha convergido a su distribución estacionaria.

En la Figura 2 se muestra el comportamiento de las cadenas para los parámetros de precipitación y temperatura mínima en la semana 1 del modelo BYM, a pesar de tener algunos cambios en el patrón, cerca de las iteraciones 15000 y 30000 por ejemplo, se puede decir que se alcanza estabilidad en la cadena, conclusión que se obtiene al ver en conjunto la figura antes mencionada y el diagnóstico de Geweke.

Figura 2: Traceplot para parámetros de Prec y Tmin en semana 1, modelo BYM. 

El modelo usado para hacer inferencia es el Besag-York-Mollie, debido a que el criterio de información de Watanabe-Akaike(WAIC) tiene un valor inferior en 6 de las 8 semanas de estudio, y en las semanas donde otro modelo ajusta mejor según el criterio del WAIC, se tiene una diferencia muy pequeña. En la Tabla 7 se muestra la medida del criterio para los 5 modelos bayesianos.

3.3.1. Resultados del modelo Besag-York-Mollie. En lo siguiente se estudia el resultado obtenido por el ajuste del modelo BYM, dado que en la sección anterior se concluye que hace un mejor ajuste de predicción según el WAIC.

Gráficos de dispersión con líneas de regresión en la Figura 3 representan las relaciones (lineales) crudas de las variables independientes con la variable dependiente para la semana 1. Dicho gráfico muestra que la tasa de incidencia de dengue (RR) está asociado positivamente con el porcentaje de pobreza en esa semana, lo mismo pasa en las semanas 2,3 y 4, mientras que para las semanas 31 a 34 cambia a negativa, con temperatura mínima se da una relación positiva en las 8 semanas, con altitud la relación es negativa en casi todas las semanas con excepción de la cuarta, donde parece no existir relación(lineal) alguna, precipitación no presenta una relación definida en la primer semana, negativa en la 32 y positiva en el resto.

El modelo BYM tiene dos efectos aleatorios, uno estructurado espacialmente, donde se incorpora la información de la matriz de vecindarios W, además de que se estima un parámetro para la varianza y otro efecto sin estructura: , independientes con distribucion normal de media nula y varianza , estimada durante el ajuste con MCMC. Este modelo es el más usado en la práctica (^{Lee, 2013}), sin embargo, se requiere estimar dos efectos aleatorios para cada dato puntual, mientras que solo la suma de estos es identificable en los datos, en la práctica uno de los efectos aleatorios es dominante sobre el otro, pero no se puede saber de antemano (^{Besag, Julian. y cols, 1991}), es por eso que el paquete CARBayes retorna la suma en las muestras del MCMC .

Figura 3: Gráficos de dispersión de variables vs riesgo relativo, semana1 

En la Tabla 8 es claro que el parámetro de varianza para el efecto no estructurado es pequeño, eso indica que el efecto es muy cercano a cero ya que la media del mismo es cero. Por otro lado, la varianza para el efecto espacial es mayor que 1 (uno) para las semanas 1, 31, 32, 33 y 34, 0.023 en la semana 4, 0.86 y 0.68 para las semanas 2 y 3 respectivamente. La fila para representa la mediana de dicho parámetro por semana, pero se debe aclarar que para cada semana existen 76 valores distintos para tal suma, ya que al municipio k-ésimo se le calcula .

Los estimadores posteriores para casos de dengue indican que las regiones de mayor incidencia en la semana 1 están localizadas al noreste y sureste de Puerto Rico, donde la región costera es la que más incidencia presenta según los resultados, mientras que los municipios con una incidencia menor se localizan en la región central del país. Por otro lado, en los residuos de la primera semana, se puede ver que no hay una estructura definida, además en un municipio se obtiene un valor residual muy grande comparado con los demás, eso puede dar a entender que existe un valor atípico en los datos para tal región, lo que puede conducir a la eliminación del mismo en el modelo, sin embargo, no se elimina pues al hacerlo se pierde mucha información, ya que no existe independencia entre los datos. En las semanas posteriores se mantiene el comportamiento de los residuos, es decir, son valores sin aparente estructura y el residuo para el municipio en cuestión sigue siendo superior con respecto a los demás, en la Figura 4 se resumen los datos de casos, las predicciones del modelo y sus respectivos residuos.

Figura 4: Casos, casos estimados y residuos de predicción en semana 1 

La variación estructural estimada después de tomar en cuenta las variables climáticas, geográfica y social, indican que los municipios con alta incidencia de dengue se localizan al sureste en las primeras 4 semanas, también al norte en las semanas 2 a 4, y al suroeste en la semana 1, para todas las semanas los municipios con mayor riesgo se encuentran cerca de la costa. La gura 5 muestra el efecto aleatorio espacialmente estructurado de las semanas 1 a 4 luego de tomar en cuenta las covariables, es decir, una representación suavizada espacialmente del riesgo residual (^{Hu, Wenbiao. y cols., 2012}).

En la Tabla 9 se muestran los parámetros de regresión para las variables climáticas, social y geográfica, aunque no se puede definir el tipo de relación entre riesgo relativo y pobreza, si es positiva o negativa, queda claro que es una variable explicativa importante para predecir el riesgo de contraer dengue en Puerto Rico.

Con excepción del resultado obtenido en la semana 32 y 34 para los parámetros de temperatura y precipitación respectivamente, el modelo de Besag sugiere que ante un aumento en lluvias o en la temperatura del país, se puede dar un aumento en la incidencia de dengue. También se puede concluir de los parámetros de altitud, que un aumento en altitud puede reflejarse en una disminución en el riesgo de enfermarse por dengue, con la excepción del estimador obtenido en la semana 4.

Figura 5: Efecto aleatorio espacial de riesgo relativo, semana 1 a semana 4 

En la Tabla 10 se presentan los intervalos de predicción para los parámetros de las covariables a un 95 % de credibilidad. Es claro en dichos intervalos que existe al menos una probabilidad de 0.95 de que el parámetro asociado a pobreza sea positivo para las semanas 1 y 2, negativo para las semanas 31 y 32, misma probabilidad para la tercera semana de obtener un parámetro positivo asociado a temperatura.

Se realizan pruebas para revisar la presencia de dependencia espacial en los residuos del modelo BYM, tanto con el índice de Geary como con el de Moran, se concluye que no hay autocorrelación espacial en dichos residuos. Nuevamente se hace uso de las funciones de spdep(^{Bivand, Roger. y cols., 2015}, 2013) en R-CRAN, geary.mc() y moran.mc(), donde se utilizan 5000 permutaciones con el objetivo de obtener mejores resultados, ya que en (^{Gelfand, Allan E. y cols., 2004}) aunque sugieren el uso de dichos índices como herramientas exploratorias de asociación espacial, también sugieren usarlos con un enfoque Monte Carlo realizando permutaciones de los datos en las regiones de estudio, y eso es precisamente lo que realizan las funciones que se usan en este apartado.

En la Tabla 11 se puede observar que para ninguna de las semanas y los índices hay evidencia suficiente para rechazar la hipótesis nula de no autocorrelación espacial en los residuos. La matriz utilizada en tales pruebas se obtiene al estandarizar por filas la matriz binaria de asociación.

Conclusiones