Modelos compartimentales.

Los modelos compartimentales (SIS, SIR, SEIR y sus variantes) son formas simplificadas de representar la evolución de una epidemia en una población grande. Estos modelos separan la población en distintos estados o “compartimentos” y tratan de explicar la forma en la que los individuos de esa población evolucionan de unos a otros. En el caso del modelo SIR (también llamado “modelo de Kermack-McKendrick”) se asume que la población (de tamaño \(N\), constante) se divide en “susceptibles”, \(S(t)\), “infectados”, \(I(t)\) y “recuperados”, \(R(t)\). El modelo asume que la tasa de disminución de la proporción de individuos susceptibles es proporcional al producto de las proporciones de susceptibles e infectados. La constante de proporcionalidad es la “tasa de contagio”. También se asume que la tasa de paso de infectado a recuperado es constante y a tal constante se la conoce como “tasa de recuperación”. Los individuos recuperados adquieren “inmunidad”, es decir, no contagian ni pueden ser contagiados. Formalmente todo esto significa que las funciones \(S,I\) y \(R\) evolucionan de aucerdo con el sistema de EDO’s \[ \mbox{(SIR)}\quad\left\{ \begin{array}{ll} \dot{S}=&-\beta \frac{SI}N,\\ \dot{I}=&\beta \frac{SI}N-\gamma I,\\ \dot{R}=&\gamma I. \end{array} \right. \] El cociente entre la tasa de contagio y la tasa de recuperación es el parámetro conocido como \(R_0\) o “número de reproducción”. Este parámetro determina el comportamiento a largo plazo de la epidemia. Esencialmente, valores de \(R_0\) menores que uno conducen a un equilibrio libre de enfermedad, mientras que si ese valor es más que uno se alcanzará, a largo plazo, un equilibrio “endémico”, es decir, con una proporción no nula de individuos infectados. Una presentación completa de estas y otras propiedades del modelo SIR, así como de otros modelos compartimentales puede consultarse en [Li, M. (2018) “An Introduction to Mathematical Modeling of Infectious Diseases”, Springer].

Cuando se trata de modelar la evolución de una epidemia que, previsiblemente, puede tener un grave impacto sobre las condiciones de vida de una sociedad organizada, en principio, el modelo SIR parece demasiado simplista. Posiblemente se pueda aceptar que la tasa de recuperación se mantenga constante en el tiempo, pero en general las sociedades actúan para reducir las tasas de contagio, introduciendo medidas de aislamiento, o a través de vacunaciones, para conseguir erradicar o, al menos, limitar el impacto de la epidemia. Hay variantes del modelo SIR que tratan de adaptarse a esta realidad. Una posibilidad es introducir una función de “mitigación”, es decir, una función que modele la forma en que la tasa de contagio va disminuyendo en el tiempo, de manera que la evolución del sistema está dado por el sistema de EDO’s \[ \mbox{(SIRM)}\quad \left\{ \begin{array}{ll} \dot{S}=&-\beta(t) \frac{SI}N,\\ \dot{I}=&\beta(t) \frac{SI}N-\gamma I,\\ \dot{R}=&\gamma I, \end{array} \right. \] donde \(\beta(t)>0\) es una función que permita modelar de forma apropiada el descenso en la tasa de contagio. Hay muchas otras posibilidades, que permiten incluir parámetros demográficos (tasas de nacimiento, muerte, inmigración o emigración) o el efecto de vacunaciones. Con vistas a un estudio rápido de la evolución de la actual epidemia de COVID19 entendemos que, en principio, las variaciones demográficas a corto plazo tendrán un efecto despreciable. Por razones que trataremos de explicar mejor más adelante esto nos ha llevado a ignorar esas posibles generalizaciones. Usar modelos SIR con vacunación puede parecer inapropiado, puesto que no hay por el momento vacunas disponibles contra el COVID-19, pero, en realidad, el paso de “susceptible” a “recuperado” (que es lo que realmente se introduce en el modelo con vacunación) puede darse de forma inadvertida en los recuentos si, como creemos, una proporción importante de casos de COVID-19 no son detectados, por manifestarse únicamente en síntomas leves, y, quizá, puede producir la inmunización (el paso a “recuperados”) de una parte de la población. Podemos modelar esta transición directa mediante el sistema \[ \mbox{(SIRV)}\quad \left\{ \begin{array}{ll} \dot{S}=&-\beta \frac{SI}N-vS,\\ \dot{I}=&\beta \frac{SI}N-\gamma I,\\ \dot{R}=&\gamma I+vS, \end{array} \right. \] La \(v\) en la etiqueta anterior hace referencia a que modelamos esta forma natural de adquirir inmunidad, aunque advertimos que el modelo SIRV no es el modelo con vacunación presente en la literatura que hemos consultado sobre el tema. Todavía hay otras variantes. En el modelo SEIR se trata de incorporar la posibilidad de que los individuos que entran en contacto con infectados pasen un periodo de latencia antes de ser trasmisores de la enfermedad. A los individuos en este estado intermedio se le denomina “expuestos”. La dinámica del modelo SEIR se describe por las siguiente EDO’s: \[ \mbox{(SEIR)}\quad \left\{ \begin{array}{ll} \dot{S}=&-\beta \frac{SI}N,\\ \dot{E}=&\beta \frac{SI}N-\alpha E,\\ \dot{I}=&\alpha E-\gamma I,\\ \dot{R}=&\gamma I, \end{array} \right. \] En el modelo SEIR el parámetro \(\alpha\) o, más bien, \(\frac 1 \alpha\) representa el tiempo medio de recuperación de un infectado.

Los modelos anteriores pueden modificarse de forma que se tengan en cuenta varios de los efectos que queremos introducir. Por ejemplo, podemos considerar un modelo SIR con mitigación e inmunización natural al mismo tiempo. A este modelo lo denominaremos SIRMV. Igualmente podemos manejar un modelo SEIRM, con una tasa \(\beta(t)\) variable en el tiempo, o un modelo SEIRV en el que sea posible la transición del estado “expuesto” a “recuperado” sin pasar por el estado “infectado”.

Análisis de regresión no lineal con modelos compartimentales.

Con todo el atractivo que puedan tener los modelos compartimentales en términos de interpretabilidad, el ajuste a series de datos reales de cualquiera de ellos requiere la búsqueda de valores apropiados de los parámetros implicados. Una posible estrategia es la calibración de los parámetros a partir de estudios independientes. Por ejemplo, podemos estimar la tasa de contagio \(\beta\) a partir de los valores disponibles en estudios de epidemias similares o de la misma epidemia en otros lugares en los que su desarrollo esté más avanzado. Esta aproximación, sin embargo, es bastante problemática. En ocasiones los parámetros incluidos en un modelo actúan de formas aproximadamente proporcionales y se produce un fenómeno parecido al de la multicolinealidad en regresión lineal. Por otro lado los valores observados de los números de casos activos son enteros y parece más lógico interpretar que se trata de versiones perturbadas de la curva \(I(t)\). Por estas razones nos proponemos asjutar modelos de regresión no lineal sobre las series de casos activos, de forma que si \(\hat{I}_t\) representa el número de casos activos en el día \(t\) entonces \[\hat{I}_t=I_\theta(t)+\varepsilon_t,\quad t=1,\ldots,T,\] donde \(I_\theta(t)\) es la función dentro del modelo correspondiente con parámetro \(\theta\) y condiciones iniciales iguales a las observadas, mientras que \(\varepsilon_1,\ldots,\varepsilon_T\) son errores centrados e independientes. Por ejemplo, en el ajuste del modelo SIR tendremos \(\theta=(\beta,\gamma)\) y la función \(I_\theta(t)\) es la (parte correspondiente a la) solución del sistema (SIR) con parámetros \(\beta\) y \(\gamma\) y condición inicial \(S_{\theta}(0)=\hat{S}_0\), \(I_{\theta}(0)=\hat{I}_0\) y \(R_{\theta}(0)=\hat{R}_0\). Asumimos conocidos estos valores iniciales, igual que la serie de casos activos \(\hat{I}_t\).

La estimación de los parámetros puede hacerse con el método de mínimos cuadrados, tal como se describe en [Bates, D.M. y Watts, D. G. (1988) Nonlinear Regression Analysis and its Applications, Wiley]. El estimador de mínimos cuadrados es el minimizador de la suma de desviaciones cuadráticas entre los valores observados \(\hat{I}_t\) y los valores ajustados del modelo, es decir, \[\hat{\theta}=\underset{\theta}{\mbox{argmin}}V(\theta),\] donde \(V(\theta)=\sum_{t=1}^T (\hat{I}_t-I_\theta(t))^2\). El valor mínimo de \(V\) es la suma de cuadrados residual, \[\mbox{RSS}=V(\hat{\theta}).\] Bajo la hipótesis de varianza constante, \(\sigma^2\), de los errores el estimador \[\hat{\sigma}^2=\frac{\mbox{RSS}}{T-p}\] (donde \(p\) es el número de parámetros del modelo) es un estimador consistente de \(\sigma^2\). A partir de \(\hat{\theta}\) y \(\hat{\sigma}^2\) se pueden calcular regiones de confianza para \(\theta\) o regiones de predicción para la respuesta.

Estudio de la evolución de la epidemia COVID-19 en algunas provincias chinas a a partir de modelos compartimentales. Comparación con el caso de Castilla y León.

Estudiamos a continuación la evolución de la epidemia COVID19 en lugares en los que se encuentra en un estadio más avanzado desde el punto de vista de los modelos compartimentales. En particular nos fijamos en los datos de algunas provincias chinas. Las siguientes gráficas muestran la evolución de los grupos de casos acumulados (azul), recuperados (verde) y muertos (rojo). En algunas de ellas se observa un corte en laa series de casos totales (al menos en Gansu y Hubei). Otras corresponden a provincias con poca población (el caso de Macau) o con pocos casos y evolución incompleta. Por ello nos centramos en las provincias con mayor número de casos con la excepción de Hubei, es decir, Anhui, Chongquing, Guangdong, Henan, Hunan, Jiangxi y Zhejiang.

En los modelos aplicados el tamaño de a población se mantiene constante en el tiempo. La expansión del COVID-19 está siendo muy rápida y por ello no parece que la hipótesis de tamaño de población constante suponga una grave desviación de la realidad. No incluimos por ello tasas de nacimiento o muerte. Además, en nuestro estudio contamos todas las salidas del estado “infectado como recuperados”, aunque una parte de ellas se correspondan con fallecimientos, pero, en todo caso, conducen a la situación de no trasmisor de la enfermedad.

Para las provincias elegidas hemos probado los modelos compartimentales descritos en la sección anterior. En los modelos que incluyen función de mitigación, ésta es tipo logístico, es decir, de la forma \[\beta(t)=\frac{\beta}{1+e^{\mu(x-x_0)}},\quad t\geq 0.\] Los parámetros \(\beta\), \(\mu\) y \(x_0\) se estiman a partir de los datos. El parámetro \(\beta\) representa la tasa de trasmisión antes de aplicar medidas de aislamiento, mientras que \(\mu\) tiene que ver con la velocidad de descenso en la tasa de transmisión y \(x_0\) es el punto central en el proceso de decrecimiento (si el decrecimiento es rápido \(x_0\) marca el “punto de cambio”). En el ajuste del modelo SIRM el valor estimado de \(x_0\) es 0 en todas las provincias en estudio, por lo que hemos incluido también en nuestro estudio el modelo SIRM0 en el que este modelo está fijado en 0, al igual que en el resto de modelos con mitigación considerados. La siguiente tabla resume la suma de cuadrados residual correspondiente al ajuste de cada modelo en cada provincia.

## SIR SIRV SIRM SIRM0 SIRMV SEIR SEIRV SEIRM
## Anhui 1534063 62965.0 86125.3 86125.3 62965.0 8299048 1540102 42028.4
## Chongqing 479192 18143.6 25108.1 25108.1 18143.6 2436145 516055 11897.8
## Guangdong 1812950 57175.1 60904.8 60904.8 57175.1 11706428 1908532 53050.0
## Henan 1632947 89088.7 117814.4 117814.4 89088.7 11233035 1672684 54556.6
## Hunan 1361477 20269.4 34238.3 34238.3 20269.4 6563321 1404415 14175.8
## Jiangxi 1133859 46182.6 63404.1 63404.1 46182.5 7391570 1191529 28227.0
## Zhejiang 2271713 71555.8 97113.2 97113.2 71555.8 10953662 2349289 73915.6

Aunque es difícil valorar la diferencia cuantitativa entre valores para los distintos modelos, es claro que el ajuste de los modelos SIR, SEIR y SEIRV es mucho pero que el de los demás modelos. Entre los modelos de tres compartimentos el modelo SIRV es el ganador en todas las provincias. El modelo SIRM0 no pierde nada frente al modelo SIRM (en el que se incluye \(x_0\) como parámetro libre). SIRM0 funciona algo peor que SIRV, pero el ajuste es razonable (luego volvemos sobre esta afirmación). El modelo SIRMV incluye como submodelos a SIRV y a SIRM. Lo que nos dice la tabla es que si el efecto de inmunización natural está incluido en el modelo entonces la función de mitigación no ayuda a mejorar el ajuste a los datos. En resumen, a la vista de la tabla los modelos de dos compartimentos que mejor describen la evolución de la epidemia son SIRV y SIRM0, por ese orden.

¿Qué significa buen ajuste o mal ajuste con estos datos? No hemos usado ningún procedimiento formal de contraste de ajuste de los modelos a estos conjuntos de datos y no creemos que sea posible extraer ninguna conclusión especialmente fuerte por esta vía con series de esta longitud. Pero una comparación gráfica de los datos de casos activos (\(\hat{I}_t\)) frente a los modelos ajustados (\(I_{\hat{\theta}} (t)\)) puede servir para tener una idea aproximada. Mostramos algunos de estos gráficos a continuación.

Estas gráficas pueden convencernos de que los modelos SIRV y SIRM0 aportan una buena descripción de la evolución de la epidemia en las distintas provincias. No se puede decir lo mismo de los modelos SIR,SEIR y SEIRV (no incluimos imágenes, pero el ajuste es en casi todos los casos muy malo). Los modelos SIRM y SIRMV tienen más complejidad (más parámetros, lo que perjudica la calidad de la estimación) y no aportan mejor ajuste.

Hasta ahora no hemos dedicado ningún comentario a los resultados obtenidos con el modelo SEIRM. En realidad es el que da menor suma de cuadrados residual en seis de las siete provincias. Pero hay un aspecto del análisis del que no hemos hablado hasta ahora. La minimizazión de la función \(V\), necesaria para calcular los estimadores \(\hat{\theta}\) no es siempre fácil de resolver. En nuestro análisis hemos usado el algoritmo implementado en la librería de . Este algoritmo combina aspectos del método de Gauss-Newton y del algoritmo de Levenberg-Marquardt y permite manejar restricciones (que imponemos para garantizar que \(\beta\), \(\gamma\), etc sean positivos). La función objetivo no es especialmente tratable y los diagnósticos de convergencia de los algoritmos han concluido que no se ha conseguido la convergencia en ninguna de las aplicaciones del modelo SEIRM a las provincias chinas. Esto nos advierte de que la función objetivo asociada a este modelo puede conducir a malas estimaciones. El objetivo de este estudio es aplicar el conocimiento adquirido en el estudio de estos datos para la predicción de la evolución de la epidemia de COVID19 en Castilla y León. Pequeñas variaciones en la estimación de los parámetros pueden conducir a variaciones brutales en las predicciones. En vista de la inestabilidad de la estimación en el modelo SEIRV preferimos basar nuestras predicciones en el ajuste de modelos SIRM y SIRV.

La pregunta que nos debemos hacer antes de usar ninguna predicción generada por esta vía es si, además de servir para la descripción de la evolución de la epidemia, los modelos tienen también alguna capacidad predictiva. Una forma de validar esto es estudiar cómo cambia la curva ajustada al aumentar el número de datos disponibles. Esto es lo que se puede ver en las siguientes animaciones.

Vemos que las curvas ajustadas cambian al avanzar al avanzar el tiempo e incorporar más datos para el ajuste. En el caso del modelo SIRV las curvas cambian de forma descontrolada en los días iniciales de la epidemia (la serie empieza el 22 de enero) y las predicciones a largo plazo no tienen ningún valor. Pero en cuanto el crecimiento empieza a ralentizarse (más o menos el 14 de febrero) la predicción se estabiliza y rápidamente se aproxima a la parte final de la serie de datos. El ajuste de modelos SIRM parece mucho más estable y, en principio, parece más fiable como base para la predicción. Por otra parte, al avanzar el tiempo las predicciones con los dos modelos tienden a ser muy parecidas. Podríamos usar al acercamiento de las dos curvas como signo de “estabilización”.

El comportamiento es parecido en otras provincias.

Los datos de Castilla y León

Analizamos la serie de datos de casos activos obtenida a partir de los informes publicados en https://analisis.datosabiertos.jcyl.es/pages/coronavirus/ La serie comienza el 1 de marzo. Las predicciones basadas en los models SIRV y SIRM evolucionan tal como se ve en las siguientes animaciones.

Como en el caso analizado antes, el comportamiento de las predicciones es bastante errático en la fase inicial, con cambios aún más bruscos en el caso de las predicciones basadas en el modelo SIRV. Las predicciones SIRM se han estabilizado en los últimos días y, un poco más lentamente, también las basadas en el modelo SIRV. Mostramos ahora las dos predicciones superpuestas para comprobar esta última afirmación.

En resumen, parece que la estabilización se está produciendo en los útimos días. De mantenerse la tendencia, a partir de ahora las predicciones serán más estable y, esperamos, más fiables. Las últimas curvas predichas son las siguientes:

Para el modelo SIRM, que es el que produce predicciones más estables, el “pico” de infectados se espera el día 7 de abril con 5582 casos activos (esta predicción se basa en los datos disponibles el 31 de marzo). A partir de ese momento se espera un descenso, con una reducción del 90% con respecto a este máximo a mediados de mayo (previsión de 574 casos activos el 15 de mayo, aunque es todavía pronto para una predicción con un plazo tan largo).

Conclusiones

El ajuste de modelos de regresión no lineal a los datos de evolución de la epidemia de COVID19 en varias provincias chinas parece mostrar que alguna variante de modelos SIR proporciona una descripción razonable de dicha evolución. Los modelos SIRV (SIR con inmunización natural) y SIRM (SIR con mitigación) son los modelos que mejor comportamiento han mostrado. Estos modelos no sirven para la predecir el momento ni la magnitud del pico de la epidemia en fases iniciales, pero al empezar la ralentización en el número de nuevos casos sí que se aproximan rápidamente a la evolución futura de la epidemia. El comportamiento de las predicciones basadas en el modelo SIRM parece más estable y, por ello, más fiable. En el caso de Castilla y León parece que hemos alcanzado ese punto de estabilización en las predicciones SIRM. Esperamos que los próximos datos confirmen esta tendencia.

Nota

Junto al enfoque anterior, se ha utilizado otro enfoque explicado brevemente en este informe que se actualiza y traslada cada 3/4 días a la Junta de CyL.