El modelo logístico y COVID-19

La modelación de crecimiento de poblaciones es una de las tareas científicas más recurrentes (y en alguna forma también lo es dentro de la ciencia de datos). Cuando nos referimos a población nos estamos refiriendo, a personas pero también puede verse como dinero, en este caso personas enfermas, etcétera.

En esta entrada abordaremos la modelación del crecimiento de COVID-19, enfermedad respiratoria causada por el virus SARS-CoV-2. Para esto ocuparemos el modelo logístico, que es el más popular y común para modelar crecimientos poblacionales. Usaremos los datos del caso italiano pero puede ser aplicable para cualquier otro conjunto de datos con las respectivas consideraciones específicas de cada elección.

Puede que haya otros modelos mucho más refinados, sin embargo la intención de esta entrada es abordar el proceso de interpretación de un modelo así como su ajuste de parámetros. No ocuparemos el mexicano debido a que todavía nos encontramos en la fase exponencial. Los paquetes de Python que ocuparemos son los siguientes:

import pandas as pd
import numpy as np
from datetime import datetime,timedelta
from sklearn.metrics import mean_squared_error
from scipy.optimize import curve_fit
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
import seaborn as sns
from pprint import pprint
%matplotlib inline

Obtención y preprocesamiento de datos.

Una de las partes más importantes de los procesos de modelación es la obtención de datos y la limpieza de los mismos. Pandas puede recuperar datos desde varias fuentes de datos, en este caso lo recuperaremos desde una url como se muestra a continuación.

url = "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv"
df = pd.read_csv(url)
df.iloc[:,[0,11]].tail()

Nos vamos a quedar con las columnas 0 (fecha)y 11 (casos_totales). La última línea nos muestra las últimas 5 líneas del dataset.

Selección_212

Ahora, si examinamos las fechas podemos ver que tienen un formato que no podemos usar así que usaremos las funciones python.datetime. Primero hacemos que nuestro dataframe sea la selección de columnas de nuestro interés (0 y 11). Después de eso pondremos las fechas en la variable date y con el formato definido "%Y-%m-%dT%H:%M:%S", que es el formato en el que vienen los dato, aplicamos la función fdt para generar toda una columna de fechas manipulables por python.

df = df.iloc[:umb,[0,11]]
formato = "%Y-%m-%dT%H:%M:%S"
date = df.data
fdt = lambda x : (datetime.strptime(x,formato) - datetime.strptime("2020-01-01T00:00:00", formato)).days
df['data'] = date.map( fdt )

Nuestro dataframe df ahora tiene dos columnas, una de fechas manipulables por Python y la otra de número de casos totales. Vamos a usar el paquete Seaborn para ver como se ven estos puntos en scatter plot.

sns.scatterplot(x=df['data'], y=df['totale_casi'])

La gráfica se ve así

Selección_213

Modelo Logístico

Hemos escuchado el término de «crecimiento exponencial» y hemos visto las gráficas que lo exhiben ya que han estado en redes sociales, medio de comunicación, vaya, en todos lados. Esto puede resultar angustiante porque efectivamente muestran el crecimiento que parece desmedido de este particular virus que tiene la característica de ser muy contagioso. Una de las características del crecimiento exponencial es que una proporción de la población, al infectarse, deja de ser susceptible y se convierte en infectada con la capacidad de transmitir el virus.

Sin embargo, este proceso no puede ser infinito, no puede crecer para siempre básicamente porque esa población susceptible va disminuyendo conforme pasa el tiempo y los que eran infectados se convierten en recuperados (o bien en decesos pero el modelo que ocuparemos hace la consideración de que la población es fija). En algún momento la capacidad del virus para infectar agentes susceptibles se detiene debido a que una buena parte de la población ya tiene memoria inmunológica que impide que el virus se propague con la misma facilidad que al inicio del proceso. La expresión que modela este comportamiento es la siguiente:

f_{a,b,c} (x) = \frac {c} {1 + \exp( -(x-b)/a)}

A esta función la llamamos modelo logístico y tiene tres parámetros. En donde 𝑥 representa la población cuyo crecimiento deseamos medir como crece, y los números 𝑎,𝑏,𝑐 son parámetros de nuestro modelo. 𝑎 mide la rapidez de la infección, 𝑏 es el día con la mayor cantidad de infecciones, 𝑐 es el número total de casos registrados al final de la infección. Decimos que c regula el final del proceso de infección ya que cuando el denominador se acerca a 𝑐 entonces 𝑓(𝑥) vale 1 que es el valor máximo que puede tomar.

Para encontrar los parámetros que mejor ajustan nuestros datos. Esto lo podemos hacer con la función scipy.optimize.curve_fit. Primero definiremos nuestra función y luego nuestras variables dependiente (x) y variable de respuesta (y).

def logistico(x, a,b,c):
    """
    Modelo logístico con argumento x (días) y con parámetros
    P=[a,b,c]
    Regresa la función f(x) = 1/(1+e(-(x-b)/a))
    """
    return c/(1 + np.exp(-(x-b)/a))

x = df.iloc[:,0]
y = df.iloc[:,1]

Para encontrar el mejor ajuste usaremos curve_fit pasándole x e y con un punto inicial para cada parámetro. Estos valores iniciales se pueden encontrar a través de otros métodos.

popt, pcov = curve_fit(logistico,x,y,p0=[2,100,20000])

En popt se guardan los valores encontrados mientras que en pcov se guarda la matriz de covarianza de los parámetros encontrados. Los parámetros encontrados son los siguientes

pprint(popt)
a, b, c = popt
array([5.22809110e+00, 8.13290036e+01, 1.20569392e+05])

La matriz de covarianza pcov nos muestra como variaron los parámetros entre sí, la diagonal es la varianza de cada parámetro. Vamos a tomar el la raíz cuadrada para calcular los errores estándar.

errores = [np.sqrt(pcov[i,i]) for i in [0,1,2]]
print(errores, sep='\t')
[0.07688318769831534, 0.22581878353538162, 2354.0038555475603]

Entonces la cantidad total de infectados se espera alrededor de c error_c

c_rango = list(map(np.round, [ popt[2] - errores[2] , popt[2] + errores[2]]))
c_rango
[118215.0, 122923.0]

El pico de la infección se espera el día b del año, es decir el día 81 del año, eso corresponde al 22 de Marzo.

pico = datetime.fromordinal(int(b))
print(b)
print(pico)
81.32900363914275
0001-03-22 00:00:00

Ahora nos gustaría saber cuándo terminará el proceso de infección y esto se puede determinar a partir cuando el acumulado de las personas infectadas iguale al parámetro 𝑐. Para esto requerimos usar la función fsolve para encontrar las raíces de 𝑓(𝑥) tomando como punto inicial el punto 𝑏

sol = int(fsolve(lambda x: logistico(x,a,b,c) - int(c),b))
sol
147

El día 147 corresponde al 27 de mayo

fin = datetime.fromordinal(sol)
fin
datetime.datetime(1, 5, 27, 0, 0)

Al momento de la publicación la cantidad de casos confirmados es de 147,577. Es importante notar que estos modelos son muy sensibles a las condiciones iniciales por lo que una variación pequeña puede resultar en cambios importantes en el resultado final.

Los datos graficados vs. la predicción del modelo lo podemos ver en la siguiente gráfica.

pred_x = list(range(max(x),sol))
plt.rcParams['figure.figsize'] = [7, 7]
preds = list(set(x).union(pred_x))
plt.rc('font', size=14)

plt.scatter(x,y,label="Real data",color="red")

plt.plot(preds, [logistico(i,popt[0],popt[1],popt[2]) for i in preds], label="Modelo logístico" )

Selección_214

En esta imagen la línea azul representa la predicción mientras que los puntos rojos representan los datos originales.

Conclusiones

Hay modelos más precisos para ajustar los datos de un proceso infeccioso. Particularmente hay modelos basados en ecuaciones diferenciales que pueden ser más precisos al tomar más consideraciones, de hecho parece que dadas las medidas de distanciamiento social la evolución de casos confirmados parece variar y se requieren modelos con parámetros distintos para diferentes etapas. Sin embargo en esta entrada lo que hicimos fue ajustar un modelo logístico, que es de los más comunes para modelar el crecimiento de una población, usando algunos parámetros que hay que determinar. En este particular caso la población es la cantidad de enfermos confirmados y los parámetros encontrados a través de resolución de ecuaciones no lineales.