En esta entrada se analiza la evolución del coronavirus en España durante las últimas   semanas, utilizando los últimos datos disponibles en Wikipedia.

El objetivo de este artículo es doble.  Además de repasar los últimos datos aparecidos, la idea es analizar el patrón de crecimiento mediante la estimación de un Generalized Growth Model (GGM) como en Chowell et al. (2017). y tratar de predecir su evolución.

La modelización matemática y estadística en epidemiología, es especialmente útil para estudiar la propagación de epidemias y puede servir para aportar evidencia sobre las trayectorias más plausibles para las próximas semanas. El ejercico de modelización que presento parece bastante sensible al último dato recogido el 9 de marzo por lo que los resultados hay que tomarlos con cautela.  Aún así, los resultados del GGM estimado con datos para el período comprendido entre el 23 de febrero y el 9 de marzo apuntan a que:

(i) el número de contagiados más plausible para final de semana es de 7.766 y

(ii) que dentro de dos semanas el número de contagios ascendería a los 88.782 casos.

Los datos

El motivo de utilizar los datos Wikipedia y no los oficiales, se debe a que éstos están más en línea con los reportados por otras webs que se están actualizando continuamente como son la de RTVE, Worldmeter o el Johns Hokpins. De hecho, como se ve en el Gráfico 1, el Ministerio de Sanidad y el European Center of Disease Control ECDC parecen estar reportando menos casos que el resto de fuentes de forma sistemática.  En cualquier caso, y con independencia de las discrepancias, en el número exacto de los casos (a las 00.00h del 10 de marzo la Wikipedia apuntaba 1.231 casos vs Ministerio/ECDC con 1.024 casos), el patrón observado en el Gráfico (1) es claramente geométrico, tal y como se mencionó en la entrada anterior aquí  a fecha de miércoles 4 de marzo.

Gráfico 1. Evolución de contagiados por Covid19 en España 24 febrero-09 marzo

El Gráfico (2) muestra el número de incidencias o casos diarios. Como puede apreciarse, en los últimos días, el número de casos se ha incrementado sustancialmente. El día 6 de marzo se añadieron 119 casos, el 7 de marzo 115 casos y el 8 de marzo 158 casos. La última actualización en el momento de escribir este artículo (a las 00.00h del 10/03/2020) implica un crecimiento durante el 9 de marzo muy abrupto, de 558  casos.

Gráfico 2. Evolución nuevos contagios 1 febrero – 09 marzo

Nota: cálculos basados en la actualización de las 00.00h del 10/03/2020, Wikipedia 

 

El Gráfico 3 a continuación, aporta información de (i) las tasas de variación diarias y de  los factores de (ii)  duplicación y  (iii) multiplicación  por 10, desde el 24 de febrero (fecha del despegue) hasta el 9 de marzo. Los datos implican que para este período, el número de contagiados se está duplicando cada 1.65 días (aprox 40 horas) y que se multiplicaría por un factor de 10 cada 5.48 días (aprox 132 horas).

Gráfico 3: Tasas de crecimiento diarias y factores de multiplicación

Nota: cálculos basados en la actualización de las 00.00h del 10/03/2020, Wikipedia 

El modelo

Aunque el cálculo de estos factores de multiplicación es informativo de la seriedad de la situación, la estimación de modelos econométricos ampliamente aceptados en el campo de la epidemiología puede permitirnos obtener conclusiones más precisas sobre la propagación del virus. En este sentido, la modelización de epidemias tiene todo un repertorio de modelos matemáticos y estadísticos para analizar su progresión tras la fecha del estallido. Algunos bastante comunes son: (i) el SIR (Susceptibles-Infectados-Recuperados) y  (ii) el SEIR (Susceptibles-Expuestos-Infectados-Recuperados), que utilizan la información de varias variables que desafortunadamente no están disponibles para el caso de España (por lo que he mirado, no existe una serie temporal consolidada de las recuperaciones).

Para analizar la evolución del número de contagiados de Covid-19 en España se emplea un Generalized Growth  Model que sólo requiere la estimación de dos parámetros clave y que se ha mostrado muy eficaz a la hora de predecir outbreaks víricos en otros contextos (Viboud et al., 2016). Una ventaja de este modelo con respecto a un modelo de crecimiento exponencial simple es que permite relajar el supuesto de crecimiento exponencial en las primeras fases mediante un parámetro de escala p. El modelo viene dado por la siguiente ecuación diferencial:

donde C'(t) describe el crecimiento de las incidencias en el período t, C(t) es el número total de contagiados en t, r es un parámetro positivo que captura la tasa de crecimiento y p es el parámetro de desaceleración del crecimiento, que se supone entre 0, 1.

Si p=0, esta ecuación describe una incidencia constante en el tiempo y el número de casos crece de forma lineal, mientras que si p=1 el modelo describe un patrón de crecimiento sub-exponencial (i.e, patrón polinómico).  La estimación de los parámetros clave del modelo r y p, utiliza los datos existentes desde el 23 de febrero al 9 de marzo y se lleva a cabo mediante NLSQ (Non-linear Least Squares).  Como se observa en el Gráfico 4, aunque el modelo es relativamente preciso a la hora de explicar el patrón temporal de los contagiados de Covid-19 tiende a infraestimar ligeramente los contagios del 1 al 4 de marzo y a sobreestimarlos en los días 7 y 8 de marzo. El error más elevado está en el último dato del 9 de marzo. Como vamos a ver, el que el cambio brusco de ayer 9 de marzo sea una nueva tendencia o bien un outlier, tiene consecuencias importantes para la predicción.

Gráfico 4: Ajuste dentro de muestra

En el Gráfico 5 a continuación, se muestra la distribución aproximada de los parámetros r y p mediante bootstrapping (asumiendo una distribución Poisson en el error)  con 200 repeticiones.  El parámetro de la tasa de crecimiento r toma el valor mediano de 0.6364 mientras que el valor del parámetro de escala p toma el valor de 0.938, lo que sugiere un comportamiento sub-exponencial.

Una de las cuestiones importantes que nos permite la estimación de los parámetros clave del GGM en la fase temprana de la epidemia es que con ellos en mano, es posible determinar el R0 del virus, es decir, el número promedio de casos nuevos que genera un caso a lo largo de un período infeccioso (Chowell et al, 2016). Sin embargo, para obtener el  R0 es preciso una estimación externa de lo que en epidemiología se denomina el intervalo de serie del virus (Tg), que se refiere al tiempo entre casos sucesivos en una cadena de transmisión. Las estimaciones existentes sobre Tg en el caso del coronavirus apuntan a que este parámetro está entre 4 y 5.2 (ver aquí y aquí) lo que indicaría un valor promedio del R0 centrado entre 2.54 y 3.3 para España, por encima de los datos de la Organización Mundial de la Salud (OMS) que estimaba que el valor del R0  se situó entre el 2 y 2.5 para el caso de China.

Gráfico 5: Distribución de los parámetros del modelo e incertidumbre

En el panel inferior, del Gráfico 5 se incorpora la incertidumbre a la predicción mostrada en el Gráfico 4 derivada del posible error de estimación de los parámetros.

Una de las cosas que podemos hacer con el modelo es utilizarlo para predecir la evolución de casos. Es importante destacar que aunque el modelo ajuste bien los datos dentro de muestra, otra cosa muy diferente es utilizarlo para predecir el futuro con certeza. Para ver si un modelo es fiable en sus pronósticos, en casi todos los campos se suelen realizar ejercicios de «backtesting», que básicamente consisten en probar la capacidad predictiva del modelo haciendo pronósticos fuera de muestra en momentos pasados y comparar los resultados con lo que finalmente sucedió. Aunque un análisis más riguroso de su fiabilidad pasaría por estimar el período en una ventana móvil e ir haciendo los pronósticos recursivamente,  el ejercicio aquí es más sencillo. Consiste en alimentar el modelo únicamente con datos del período comprendido entre el 23 de febrero y el  4 marzo y ver, si las predicciones para los días del 5, 6, 7 y  8 de marzo no se desvían mucho de lo que ha sucedido. El Gráfico 6 muestra que el resultado de este ejercicio de backtest: los intervalos de confianza de la predicción a 4 días vista cubren en todo momento la trayectoria observada. Por tanto, el miércoles pasado ya se podía prever que para el domingo 8 como poco habría 600 casos y como máximo 1200 casos.

Gráfico 6: Backtesting del GGM,  4 marzo – 8 marzo

Finalmente, el Gráfico 7 muestra los resultados de utilizar el modelo con los parámetros estimados r=0.6364 y p=0.938, teniendo en cuenta su incertidumbre (outputs de los Gráficos 4 y 5) y utilizando toda la información disponible desde el 23 de febrero al 9 de marzo para pronosticar la trayectoria del 10 al 24 de marzo. Por otro lado, en la Tabla 1 se reporta la predicción mediana o central por días, así como el intervalo de confianza del 95%. Los resultados apuntan a que con un 95% de probabilidad los contagios se situarán entre los 6400-9991 siendo el escenario más probable el de los 7.766 contagios.   Para dentro de dos semanas el escenario central sería de 88.782 casos y el intervalo de [56.300-155.776] casos. Como se observa en el ensanchamiento de las bandas de confianza del Gráfico 7,  la incertidumbre es muy elevada y  crece considerablemente conforme nos alejamos del origen del pronóstico.

Gráfico 7: Predicciones modelo

Tabla 1: Proyecciones del número de contagios en España

En cualquier ejercicio de modelización y/o pronóstico un punto importante es el de la robustez. En este sentido, es importante mencionar que el dato de ayer del 9 de marzo es bastante clave en la generación de este resultado tan explosivo, y que los resultados presentados NO son robustos. De hecho, el mismo modelo estimado y calibrado únicamente hasta el 8 de marzo (adjunto las trayectorias pronosticadas dentro y fuera de muestra abajo) no implicaba una evolución tan explosiva. A una semana y dos semanas vista, las horquillas que cubren el 95% de probabilidad serían de [5-9 mil] y de [15-34 mil] contagiados respectivamente.

El que el dato del 9 de marzo finalmente sea un outlier o un cambio drástico de tendencia, va a ser clave para sacar conclusiones y lo sabremos pronto. En cualquier caso, tanto si nos quedamos con el output del modelo alimentado únicamente con información hasta el 8 de marzo  (mínimo 4 mil casos adicionales) como si nos fiamos de lo que dice el modelo alimentado con información más reciente (mínimo 5200 casos adicionales),  los resultados de este análisis sugieren que las medidas de control de la epidemia en España tendrán que ser muy enérgicas si no queremos que el coronavirus se extienda masivamente por la población y el territorio.

En el caso de Italia las medidas de cierre de colegios en todo el país se tomaron hace ya días 5 días con 3 mil contagios y aún con ello, el número ha crecido a 9 mil infectados. Está por ver su efectividad pero ayer se tomaron medidas drásticas de cuarentena y restricción a la movilidad en todo el territorio italiano hasta el 4 de abril. Aunque sean duras, parecen correctas puesto que existe evidencia que apunta a que las restricciones a la movilidad implementadas en China sirvieron para evitar los escenarios de crecimiento exponencial. Dejo los enlaces aquí y aquí.

—————————————————————————————————————————————
*** Actualización a las 15.00h del 10/03/2020: Con 1592 casos se confirma que el dato del 9 de marzo no era un outlier ya que el crecimiento observado parece ajustarse al modelo usando toda la información muestral. De hecho, puede que a final del día el dato sea tan elevado que caíga en el 2.5% de probabilidad que tenían las trayectorias de contagios por encima de los 1748 casos.

*** Actualización a las 00.00h del 11/03/2020: Con 1695 casos se confirma que el dato del 9 de marzo no era un outlier ya que el crecimiento observado parece ajustarse al modelo usando toda la información muestral.  El resultado del pronóstico a 1 día finalmente cae dentro del intervalo del 95%.