En este artículo analizamos los datos oficiales de la evolución de la pandemia de COVID-19 en México mediante un modelo de mapeo de enfermedades de Poisson. La tasa de nuevos casos, relativa a la población en riesgo, se descompone en efectos factoriales definidos por variables categóricas como sexo, grupo de edad, entidad federativa y semana epidemiológica. Estos efectos factoriales nos permiten intepretar distintas dimensiones de la pandemia en México. La estimación de los efectos se realiza de manera bayesiana, cuantificando así toda la incertidumbre a través de distribuciones de probabilidad.
En Diciembre de 2019, la ciudad de Wuhan, de la provincia Hubei en China, se convirtió en el centro de un brote de neumonía de causa desconocida, lo cual atrajo la atención internacional [Wan+20]. Un grupo de pacientes se presentó a diferentes hospitales con diagnósticos de neumonía de etiología no conocida. La mayoría de estos pacientes fueron vinculados epidemiológicamente a un mercado mayorista de pescados, mariscos y animales vivos en la provincia de Hubei [Bog+20]. Autoridades chinas de salud hicieron una investigación para caracterizar y controlar la enfermedad.
En pocos días los contagios aumentaron exponencialmente, no sólo en la parte continental de China sino también en diferentes países. El agente causal fue identificado; un nuevo coronavirus (2019-nCoV) posteriormente clasificado como SARS-CoV2 causante de la enfermedad denominada como COVID-19. El 11 de marzo del 2020 la Organización Mundial de Salud declara esta enfermedad como una pandemia [KH20].
De acuerdo con [Suá+20], el primer caso confirmado de COVID-19 en México fue el 27 de febrero de 2020. Casi cuatro semanas después, el 23 de marzo, la Secretaría de Salud en México reconoció al COVID-19 como una enfermedad de prioridad nacional. Desde esta fecha y hasta el 30 de marzo, cuando la situación fue declarada emergencia de salud nacional, el número de casos había alcanzado ya 2,694. La estrategia nacional consistió en comunicación masiva, campañas de prevención y la suspensión inmediata de actividades económicas no esenciales desde el 30 de marzo y hasta el 1 de junio.
Desde el comienzo de la pandemia en nuestro país, la Secretaría de Salud en México produce un reporte diario de casos confirmados de COVID-19. La confirmación se realiza mediante la aplicación de una prueba de reacción en cadena de la polimerasa (PCR, por sus siglas en inglés) [e.g. WWB20], una prueba de antígeno y, como en el caso de las defunciones, mediante un comité dictaminador para casos en los que no se pudo tomar una muestra o bien la muestra fue inválida. Dicho reporte diario contiene una lista de los casos confirmados acumulados hasta la fecha de publicación, en el cual se incluyen la entidad federativa, el sexo, la edad y la fecha de inicio de síntomas.
El objetivo de este trabajo es analizar los datos de casos confirmados de COVID-19 en México y descomponerlos en efectos factoriales. Esto se realiza mediante un modelo de mapeo de enfermedades con base en una regresión de Poisson. Se definen la tasa del total de nuevos casos relativos a la población en riesgo, y se descompone en efectos factoriales definidos por variables categóricas como sexo, grupo de edad, entidad federativa y semana epidemiológica en la cual iniciaron los síntomas.
La estructura del resto del artículo es la siguiente: en la Sección 2 se describe el modelo de regresión de Poisson y la descomposición en efectos factoriales. En la sección 3 se presentan las distribuciones iniciales necesarias para la inferencia bayesiana. En la Sección 4 se presentan los resultados. Finalmente concluimos el trabajo en la Sección 5.
Antes de proceder, introducimos la notación. $\po(\mu)$ denota una distribución de Poisson con media $\mu$; $\no(\mu,\tau)$ denota una distribución normal con media $\mu$ y precisión $\tau$ (i.e. varianza $1/\tau$); $\ga(\alpha,\beta)$ denota una distribución gamma con media $\alpha/\beta$.
Los modelos de mapeo de enfermedades tienen por objetivo el describir y cuantificar la variación espacial en el riesgo de una enfermedad, incluyendo la identificación de patrones de posible asociación entre vecinos [Law18]. Por lo cual, como resultado del estudio, es posible identificar áreas de alto o bajo riesgo para contribuir al entendimiento de la enfermedad.
En nuestro estudio de los casos oficiales de COVID-19 en México usaremos cuatro variables categóricas: sexo con categorías de $1=$ Hombre y $2=$ Mujer; grupo de edad con categorías en intervalos de años $1=[0,20]$, $2=(20,40]$, $3=(40,60]$ y $4=(60,\infty)$; entidad federativa con categorías dadas por cada uno de los 32 estados de México según la ordenación oficial del Instituto Nacional de Estadística y Geografía (INEGI); y la semana epidemiológica de inicio de síntomas donde $1=$ Semana 10, hasta $38=$ Semana 47.
Sean $Y_{i,j,k,l}$ el número de casos nuevos de COVID-19 y $X_{i,j,k}$ la población en riesgo (en cientos de miles) de sexo $i$, grupo de edad $j$, entidad federativa $k$ y semana epidemiológica $l$, para $i=1,2$, $j=1,\ldots,4$, $k=1,\ldots,32$ y $l=1,\ldots,38$. Nótese que la población es considerada constante a lo largo de las semanas epidemiológicas *La información poblacional es parte de las proyecciones poblacionales del Consejo Nacional de Población (CONAPO) http://www.conapo.gob.mx/es/CONAPO/.
Consideramos el siguiente modelo de regresión de Poisson [e.g. NA14]
\begin{equation} \label{eq:mod} Y_{i,j,k,l}\sim\po\left(X_{i,j,k}\lambda_{i,j,k,l}\right), \end{equation}donde $\lambda_{i,j,k,l}$ es la tasa de incidencia COVID-19 por cada cien mil habitantes de sexo $i$, en el grupo de edad $j$, en la entidad federativa $k$ y en la semana epidemiológica $l$. Esta tasa, en escala logarítmica, se define a través de los efectos principales asociados a los niveles de cada una de las variables categóricas o factores:
\begin{equation} \label{eq:tasa} \log\left(\lambda_{i,j,k,l}\right)=\alpha+\beta_i+\gamma_j+\delta_k+\epsilon_l, \end{equation}donde $\alpha$ corresponde al nivel medio global; $\beta_i$ es el efecto del sexo $i$, $i=1,2$; $\gamma_j$ es el efecto del grupo de edad $j$, $j=1,\ldots,4$; $\delta_k$ es el efecto de la entidad federativa $k$, $k=1,\ldots,32$; y $\epsilon_l$ es el efecto de la semana epidemiológica $l$, $l=1,\ldots,38$.
Para que el efecto medio global sea estimable, es necesario imponer una serie de restricciones a los demás parámetros. La restricción común es hacer el efecto de la primera categoría de todas las variables categóricas igual a cero, i.e., $\beta_1=\gamma_1=\delta_1=\epsilon_1=0$. Sin embargo, esta restricción usualmente complica la interpretación de los demás parámetros y confunde las categorías base, primer nivel de cada variable categórica, con el efecto global. Una restricción alternativa que facilita la interpretación de los parámetros es la siguiente:
\begin{equation} \label{eq:restriccion} \sum_{i=1}^2\beta_i=\sum_{j=1}^4\gamma_j=\sum_{k=1}^{32}\delta_k=\sum_{l=1}^{38}\epsilon_l=0. \end{equation}La restricción \eqref{eq:restriccion} permite que $\alpha$ se interprete como el nivel medio global y cada uno de los demás coeficientes como cambios respecto al nivel medio. En este trabajo usaremos esta última restricción.
La descomposición factorial de la tasa \eqref{eq:tasa} permite generar tasas parciales. Por ejemplo, $\exp\{\alpha\}$ es la tasa de incidencia media por COVID-19 en México, $\exp\{\alpha+\beta_i\}$ es la tasa de incidencia por COVID-19 para el sexo $i$, $\exp\{\alpha+\gamma_j\}$ es la tasa de incidencia por COVID-19 para el grupo de edad $j$, $\exp\{\alpha+\delta_k\}$ es la tasa de indicencia por COVID-19 en la entidad federativa $k$, y $\exp\{\alpha+\epsilon_l\}$ es la tasa de incidencia por COVID-19 en la semana epidemiológica $l$.
Sea $\bY=\{Y_{i,j,k,l}\}$ una colección de variables aleatorias provenientes del modelo \eqref{eq:mod} y sea $\btheta=\{\alpha,\beta_i,\gamma_j,\delta_k,\epsilon_l\}$ el vector de parámetros desconocidos del modelo, para $i=1,2$, $j=1,\ldots,4$, $k=1,\ldots,32$ y $l=1,\ldots,38$. Entonces toda la información disponible de la muestra se resume en la función de densidad conjunta (función de verosimilitud) dada por $$f(\by\mid\btheta)=\prod_{i,j,k,l}\po\left(y_{i,j,k,l}\mid\exp\{\alpha+\beta_i+\gamma_j+\delta_k+\epsilon_l\}X_{i,j,k}\right).$$ Es posible obtener los estimadores máximo verosímiles de $\btheta$ maximizando la función de verosimilitud anterior sujeto a las restricciones \eqref{eq:restriccion}. Alternativamente podemos usar un enfoque bayesiano para realizar la inferencia, para lo cual es necesario describir nuestro conocimiento inicial sobre $\btheta$ a través de distribuciones de probabilidad. Las distribuciones iniciales elegidas son: $\alpha\sim\no(0,0.001)$; $\beta_i\mid\tau_\beta\sim\no(0,\tau_\beta)$ para $i=1,2$ y $\tau_\beta\sim\ga(0.1,0.1)$; $\gamma_j\mid\tau_\gamma\sim\no(0,\tau_\gamma)$ para $j=1,\ldots,4$ y $\tau_\gamma\sim\ga(0.1,0.1)$; $\delta_k\mid\tau_\delta\sim\no(0,\tau_\delta)$ para $k=1,\ldots,32$ y $\tau_\delta\sim\ga(0.1,0.1)$; $\epsilon_1\mid\tau_\epsilon\sim\no(0,\tau_\epsilon)$, $\epsilon_l\mid\tau_\epsilon\sim\no(\epsilon_{l-1},\tau_\epsilon)$ para $l=2,\ldots,38$ y $\tau_\epsilon\sim\ga(0.1,0.1)$. Suponiendo independencia entre parámetros de distinto tipo tenemos $$f(\btheta)=\no(\alpha\mid 0,0.001)\prod_{i=1}^2\no(\beta_i\mid 0,\tau_\beta)\prod_{j=1}^4\no(\gamma_j\mid 0,\tau_\gamma)\prod_{k=1}^{32}\no(\delta_k\mid 0,\tau_\delta)$$ $$\hspace{-4.8cm}\times\no(\epsilon_1\mid 0,\tau_\epsilon)\prod_{l=2}^{38}\no(\epsilon_l\mid\epsilon_{l-1},\tau_\epsilon)$$ $$\hspace{1.3cm}\times\ga(\tau_\beta\mid 0.1,0.1)\ga(\tau_\gamma\mid 0.1,0.1)\ga(\tau_\delta\mid 0.1,0.1)\ga(\tau_\epsilon\mid 0.1,0.1)$$
La distribución inicial para $\alpha$ es del tipo vaga o difusa, en el sentido de que representa un conocimiento inicial vago que tiene poco impacto en la inferencia posterior. La distribución inicial para los parámetros $\{\beta_i\}$, $\{\gamma_j\}$ y $\{\delta_k\}$ son del tipo jerárquicas, es decir, imponen cierto grado de dependencia entre parámetros del mismo tipo, e.g. $\Cov(\beta_i,\beta_{i'}) \mayorque 0$ para $i\neq i'$, lo que implica que la información proveniente de la muestra se compartirá entre dichos parámetros, reduciendo así la incertidumbre. Por último, la distribución inicial de $\{\epsilon_l\}$, que corresponde a los efectos de la semana epidemiológica, es del tipo dinámico que impone una dependencia del tipo autorregresiva de orden uno.
Tanto la información proveniente de la muestra, resumida en la función de verosimilitud $f(\by\mid\btheta)$, como la información inicial $f(\btheta)$, se combinan a través del Teorema de Bayes produciendo un conocimiento posterior $f(\btheta\mid\by)$, es decir, $f(\btheta\mid\by)=f(\by\mid\btheta)f(\btheta)/f(\by)$, en la cual $f(\by)=\int f(\by\mid\btheta)f(\btheta)\mathrm{d}\btheta$.
La obtención de la distribución posterior de manera analítica no es posible. Sin embargo, existen actualmente métodos numéricos basados en simulaciones, llamados métodos Monte Carlo via Cadenas de Markov (MCMC) [e.g. GL06], los cuales permiten obtener muestras de $f(\btheta\mid\by)$ de manera iterativa. Los métodos MCMC para modelos paramétricos estándar, como el modelo presentado aquí, se pueden implementar en paquetes de cómputo estadístico como R [Tea20] mediante el uso de las librerías Openbugs o Jags. Las inferencias que se mostrarán más adelante fueron obtenidas de esta manera.
La fuente de los datos es la página de datos abiertos del gobierno de México*https://www.gob.mx/salud/documentos/datos-abiertos que provee el Sistema Nacional de Vigilancia Epidemiológica (SINAVE). Para fines de este estudio, nos concentraremos en casos confirmados de COVID-19 con fecha de inicio de síntomas a partir del 1 de marzo, que corresponde al inicio de la semana epidemiológica 10, y hasta el 21 de noviembre, que corresponde al fin de la semana epidemiológica 47.
En el período de estudio se cuenta con $1,109,274$ casos confirmados, de los cuales 51% son hombres y 49% son mujeres; 5% son del grupo de edad entre 0 y 20 años, 40% del grupo de edad entre 20 y 40 años, 38% del grupo de edad entre 40 y 60 años y 18% corresponden al grupo de edad de mayor a 60 años; las entidades federativas con mayor número de casos son Ciudad de México con 22%, Estado de México con 7%, y Nuevo León con 6%. El estado que menos casos tiene es Campeche con 0.6% del total de casos.
Implementamos el modelo descrito en las secciones anteriores y estimamos los parámetros mediante un enfoque bayesiano. Para ello implementamos un MCMC con 5,000 iteraciones, descartando las primeras 500 como período de calentamiento y quedándonos con una simulación cada dos para disminuir la dependencia entre simulaciones y con esta información procedemos a estimar los parámetros. Para obtener los estimadores puntuales usamos una función de pérdida valor absoluto, lo que nos indica que el estimador óptimo es la mediana de la distribución posterior [e.g. BS00]. Los intervalos de credibilidad al 95% son obtenidos mediante los cuantiles 2.5% y 97.5% de la distribución posterior. Este intervalo resultad bastante angosto, lo cual se debe a que la información con la que se cuenta es muy vasta, lo que produce estimaciones muy precisas.
El estimador puntual del efecto global es $\widehat{\alpha}=2.42$, lo que implica que la tasa de incidencia media por COVID-19 en México es $e^{\widehat{\alpha}}=11.25$ casos por semana por cada cien mil habitantes, con un intervalo de credibilidad del 95% entre $11.16$ y $11.32$ casos. Con respecto al sexo, los estimadores puntuales de los efectos son $\widehat{\beta_1}=0.0589$ y $\widehat{\beta_2}=-0.0589$. Nótese que estos estimadores suman cero debido a la restricción \eqref{eq:restriccion}. Esto implica que la incidencia en hombres es $e^{\widehat{\beta_1}}-1\approx6%$ mayor que la incidencia media nacional y en las mujeres es $6%$ menor. En otras palabras, hay $e^{\widehat{\alpha}+\widehat{\beta_1}}=11.9$ hombres por semana por cada cien mil hombres y $10.6$ mujeres por semana por cada cien mil mujeres.
Con respecto al grupo de edad, los estimadores puntuales son $\widehat{\gamma_1}=-1.666$, $\widehat{\gamma_2}=0.4157$, $\widehat{\gamma_3}=0.6701$ y $\widehat{\gamma_4}=0.5803$. Nuevamente estos estimadores suman cero debido a la restricción \eqref{eq:restriccion}, pero sólo uno de ellos es negativo y los demás positivos. Esto implica que la incidencia en el grupo de edad de 0 a 20 años ($e^{\widehat{\gamma_1}}-1$) es $81%$ menor a la incidencia media nacional, lo que confirma que los niños y adolescentes casi no se contagian. En cambio, en los demás grupos de edad se tiene que en el grupo entre 20 y 40 años la incidencia es $52%$ mayor a la incidencia media nacional, en el grupo entre 40 y 60 años la incidencia es $95%$ mayor a la incidencia media nacional y en el grupo de más de 60 años la incidencia es $79%$ mayor a la incidencia media nacional. Esto confirma que los grupos de mayor incidencia son los mayores a 40 años. Al convertir estos números a tasas $e^{\widehat{\alpha}+\widehat{\gamma_j}}$, obtenemos $2.13$, $17.04$, $21.97$ y $20.09$ casos por cada cien mil habitantes en los grupos de edad $[0,20]$, $(20,40]$, $(40,60]$ y $(60,\infty)$ respectivamente.
En cuanto a la entidad federativa, Coahuila es el estado con la menor incidencia, con una tasa de $2.02$ casos por cada cien mil habitantes ($82%$ menor a la incidencia media nacional), seguido de Morelos y Estado de México con tasas de $5.03$ y $6.02$ casos por cada cien mil habitantes ($55%$ y $46%$ menores a la incidencia media nacional), respectivamente. Por otro lado, los estados con mayor incidencia son la Ciudad de México, Baja California Sur y Tabasco, con tasas de $33.82$, $24.75$ y $20.77$ casos por cada cien mil habitantes, lo que equivale a aumentos de $200%$, $120%$ y $85%$, respectivamente, relativos a la incidencia media nacional. En la Figura 1 se presentan las tasas estimadas $e^{\widehat{\alpha}+\widehat{\delta_k}}$ para las 32 entidades federativas. Además de lo ya mencionado, podemos observar estados como Chiapas y Tabasco con tasas de incidencia altas y una región central (verde oscuro) alrededor de la Ciudad de México con tasas de incidencia por debajo de la media nacional.
Finalmente, el efecto temporal dado por la semana epidemiológica, $\epsilon_l$, se aprecia mejor mediante una gráfica. En la Figura 2 se presentan las tasas semanales estimadas $e^{\widehat{\alpha}+\widehat{\epsilon_l}}$ para las semanas de la 10 a la 47 ($l=1,\ldots,38$). Adicionalmente se incluye una línea horizontal que representa la tasa media nacional. La tendencia fue creciente desde la semana 10 (que comenzó el 1 de marzo) hasta la semana 29 (que terminó el 18 de julio), llegando al pico máximo dicha semana 29 con casi 30 casos por cada cien mil habitantes. Posteriormente las tasas decrecieron hasta la semana 38 (que terminó el 19 de septiembre) y a partir de esta semana las tasas volvieron a crecer hasta la semana 47 (que terminó el 21 de noviembre y es el fin del periodo de observación), pero sin alcanzar de nuevo el pico máximo que se observó en la semana 29. Notamos también que a partir de la semana 20 (que comenzó el 10 de mayo) las tasas superaron la tasa media nacional.
La gráfica de la Figura 2 se pude interpretar como la curva temporal de evolución de la pandemia en México. Cuando esta curva presenta un decremento (e.g. después de la semana 29) se dice que la curva de casos acumulados se aplanó
. De hecho el pico se alcanzó dos semanas después de que el semáforo epidemiológico cambió de rojo a naranja en la Ciudad de México.
Los modelos de regresión de Poisson, como el utilizado en este trabajo, son útiles para el modelado de datos de conteos. Sin embargo, si los datos presentan valores en otro tipo de escalas, e.g. continua no acotada, continua acotada, conteos acotados, o categórica, es posible utilizar algún otro modelo de la clase de modelos lineales generalizados [e.g. NA14]. La descomposición en factores de los datos oficales permitió identificar otras dimensiones de la pandemia en México, como la incidencia por sexo, grupo de edad, entidad federativa y semana epidemiológica. Esta descomposición es relevante porque puede ser muy útil para priorizar estrategias de prevención o mitigación de incidencia futura en México.
El modelo utilizado, representado por las ecuaciones \eqref{eq:mod} y \eqref{eq:tasa}, considera únicamente los efectos principales, sin embargo, es posible considerar extensiones al modelo que representen más fielmente a los datos al incluir interacciones entre los efectos principales, ya sean de orden 2, e.g. sexo y edad, o de orden mayor, e.g. sexo, edad y entidad federativa. El incluir interacciones hace que la interpretación sea más específica, ya que es posible, por ejemplo, cuantificar el efecto de mujer en grupo de edad de más de 60 años, que puede ser distinto al efecto de hombre en el mismo grupo de edad. No obstante, los modelos con interacciones hacen más lento y engorroso el proceso de estimación al aumentar drásticamente el número de parámetros.
La inferencia bayesiana de este tipo de modelos es bastante accesible gracias a la creación de paquetes computacionales tipo BUGS (Bayesian Inference Using Gibbs Sampler) que pueden ser fácilmente implementados a través de R. Tanto los datos utilizados para este análisis como el código BUGS pueden ser solicitados a los autores.
Los autores agradecen el apoyo de la Asociación Mexicana de Cultura, A.C.
