Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

Aulas do Curso de Modelagem matemática IV da FGV-EMAp

Views: 2412
License: GPL3
Image: default
Kernel: SageMath 9.3

Combinando Modelos e dados: COVID-19

Ao se modelar uma epidemia real, a concordância do modelo e os dados observados é de extrema importância.

Neste Notebook vamos estudar o Modelo SEIAHR proposto para a COVID-19 por Coelho et al. Neste modelo, temos os seguintes compartimentos: diagrama

Matemáticamente: dSdt=−λ[(1−χ)S],dEdt=λ[(1−χ)S]−αE,dIdt=(1−p)αE−δI−ϕI,dAdt=pαE−γA,dHdt=ϕI−(ρ+μ)H,dRdt=δI+ρH+γA,\begin{align} \frac{dS}{dt}&=-\lambda [(1-\chi) S],\\ \frac{dE}{dt}&= \lambda [(1-\chi) S]-\alpha E,\\ \frac{dI}{dt}&= (1-p)\alpha E - \delta I -\phi I,\\ \frac{dA}{dt}&= p\alpha E - \gamma A,\\ \frac{dH}{dt}&= \phi I -(\rho+\mu) H,\\ \frac{dR}{dt}&= \delta I + \rho H+\gamma A, \end{align}

onde λ=β(I+A)\lambda=\beta (I+A).

Para facilitar nosso trabalho com dados dentro do ambiente Sage, precisamos instalar o Pandas para isso em um terminal use o seguinte comando: 

sage -sh

este comando iniciará uma subsessão do sage onde você poderá instalar o pandas com o comando usual:

pip install pandas
pip install parameter-sherpa
exit

o comando exit após a instalação visa sair da subsessão do Sage.

import numpy as np
import pandas as pd
# %display typeset

def model(t, y, params):
    S, E, I, A, H, R, C, D = y
    chi, phi, beta, rho, delta, gamma, alpha, mu, p, q, r = params
    lamb = beta * (I + A)
    # Turns on Quarantine on day q and off on day q+r
    chi *= ((1 + np.tanh(t - q)) / 2) * ((1 - np.tanh(t - (q + r))) / 2)
    return [
        -lamb * ((1 - chi) * S),  # dS/dt
        lamb * ((1 - chi) * S) - alpha * E,  # dE/dt
        (1 - p) * alpha * E - delta * I - phi * I,  # dI/dt
        p * alpha * E - gamma * A,
        phi * I - (rho + mu) * H,  # dH/dt
        delta * I + rho * H + gamma * A,  # dR/dt
        phi * I,  # (1-p)*alpha*E+ p*alpha*E # Hospit. acumuladas
        mu * H  # Morte acumuladas
    ]           

chi = .3
phi = 0.012413633926076584
beta = 0.27272459855759813
rho = 0.2190519831830368
delta = 0.04168480042146949
gamma = 0.04
alpha =  0.3413355572047603
mu = 0.02359234606623134
p = 0.7693029079871165
q = 50
r = 55

T = ode_solver()
T.function = model
T.algorithm='rk8pd'
inits = [.99, 0, 1e-4, 0, 0, 0, 0, 0]
tspan = [0,200]
T.ode_solve(tspan, inits, num_points=200, params=[chi,phi,beta,rho,delta,gamma,alpha,mu,p,q,r])

def get_sim_array(sol):
    sim = np.array([y for t,y in sol])
    return sim
get_sim_array(T.solution).shape
(201, 8)
popRJ = 6.32e6
def plot_sol(sol):
    sim = get_sim_array(sol)*popRJ
    P = list_plot(sim[:,0],legend_label='S')
    colors = ['blue','red','pink','green','yellow','orange','black','purple']
    for i,var in enumerate(['E','I','A','H','R','C','D']):
        P += list_plot(sim[:,i+1],color=colors[i+1],legend_label=var)
   
    show(P)
plot_sol(T.solution)
Image in a Jupyter notebook
sims = get_sim_array(T.solution)
sims[-1,-2]*popRJ
328307.38107624033

Obtendo os dados

Os dados usados aqui foram obtidos do site Brasil.io. Apenas os dados do estado do Rio de Janeiro foram extraídos e salvos em um CSV. 

def load_data(state):
    df = pd.read_csv(f'dados_{state}.csv')
    df['data'] = pd.to_datetime(df.data)
#     df.set_index('data', inplace=True)
    return df

dfRJ = load_data('RJ')
ld = len(dfRJ)
html(dfRJ.tail().to_html())
data date last_available_confirmed last_available_deaths incidencia_casos incidencia_morte ew
170 2020-08-22 2020-08-22 210464 15267 65.0 3428.0 34
171 2020-08-23 2020-08-23 210948 15292 25.0 484.0 35
172 2020-08-24 2020-08-24 211360 15392 100.0 412.0 35
173 2020-08-25 2020-08-25 214003 15560 168.0 2643.0 35
174 2020-08-26 2020-08-26 214003 15560 0.0 0.0 35 
subnot=1
dfRJ.set_index('data')[['last_available_confirmed','last_available_deaths']].plot();
Image in a Jupyter notebook
dfRJ['last_available_deaths'].plot()
<AxesSubplot:>
Image in a Jupyter notebook

Ajustando o Modelo aos dados

Há muitas formas de se ajustar um modelo dinâmico aos dados disponíveis. Vamos começar por meio de otimização, buscando os valores de parâmetros que minimizam o desvio entre o modelo e os dados. O ajuste será feito simultaneamente para casos e mortes acumuladas. Abaixo vamos importar a biblioteca Sherpa:

!pip install parameter-sherpa
Collecting parameter-sherpa Downloading parameter-sherpa-1.0.6.tar.gz (513 kB) |████████████████████████████████| 513 kB 4.5 MB/s eta 0:00:01 Requirement already satisfied: pandas>=0.20.3 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from parameter-sherpa) (1.3.1) Collecting pymongo>=3.5.1 Downloading pymongo-3.12.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (531 kB) |████████████████████████████████| 531 kB 7.4 MB/s eta 0:00:01 Requirement already satisfied: numpy>=1.8.2 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from parameter-sherpa) (1.19.5) Requirement already satisfied: scipy>=1.0.0 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from parameter-sherpa) (1.5.4) Requirement already satisfied: scikit-learn>=0.19.1 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from parameter-sherpa) (0.24.2) Collecting flask>=0.12.2 Downloading Flask-2.0.1-py3-none-any.whl (94 kB) |████████████████████████████████| 94 kB 2.8 MB/s eta 0:00:01 Collecting GPyOpt>=1.2.5 Downloading GPyOpt-1.2.6.tar.gz (56 kB) |████████████████████████████████| 56 kB 4.8 MB/s eta 0:00:01 Collecting enum34 Downloading enum34-1.1.10-py3-none-any.whl (11 kB) Requirement already satisfied: matplotlib in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from parameter-sherpa) (3.3.4) Collecting Jinja2>=3.0 Using cached Jinja2-3.0.1-py3-none-any.whl (133 kB) Collecting itsdangerous>=2.0 Downloading itsdangerous-2.0.1-py3-none-any.whl (18 kB) Requirement already satisfied: click>=7.1.2 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from flask>=0.12.2->parameter-sherpa) (8.0.1) Collecting Werkzeug>=2.0 Downloading Werkzeug-2.0.1-py3-none-any.whl (288 kB) |████████████████████████████████| 288 kB 12.2 MB/s eta 0:00:01 Collecting GPy>=1.8 Downloading GPy-1.10.0.tar.gz (959 kB) |████████████████████████████████| 959 kB 12.8 MB/s eta 0:00:01 Requirement already satisfied: six in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from GPy>=1.8->GPyOpt>=1.2.5->parameter-sherpa) (1.15.0) Collecting paramz>=0.9.0 Downloading paramz-0.9.5.tar.gz (71 kB) |████████████████████████████████| 71 kB 7.1 MB/s eta 0:00:01 Requirement already satisfied: cython>=0.29 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from GPy>=1.8->GPyOpt>=1.2.5->parameter-sherpa) (0.29.21) Collecting MarkupSafe>=2.0 Downloading MarkupSafe-2.0.1-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl (30 kB) Requirement already satisfied: python-dateutil>=2.7.3 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from pandas>=0.20.3->parameter-sherpa) (2.8.0) Requirement already satisfied: pytz>=2017.3 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from pandas>=0.20.3->parameter-sherpa) (2020.4) Requirement already satisfied: decorator>=4.0.10 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from paramz>=0.9.0->GPy>=1.8->GPyOpt>=1.2.5->parameter-sherpa) (4.4.2) Requirement already satisfied: threadpoolctl>=2.0.0 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from scikit-learn>=0.19.1->parameter-sherpa) (2.2.0) Requirement already satisfied: joblib>=0.11 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from scikit-learn>=0.19.1->parameter-sherpa) (1.0.1) Requirement already satisfied: pillow>=6.2.0 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from matplotlib->parameter-sherpa) (8.1.2) Requirement already satisfied: kiwisolver>=1.0.1 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from matplotlib->parameter-sherpa) (1.0.1) Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from matplotlib->parameter-sherpa) (2.4.7) Requirement already satisfied: cycler>=0.10 in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from matplotlib->parameter-sherpa) (0.10.0) Requirement already satisfied: setuptools in /home/fccoelho/Downloads/SageMath/local/lib/python3.9/site-packages (from kiwisolver>=1.0.1->matplotlib->parameter-sherpa) (51.1.1) Building wheels for collected packages: parameter-sherpa, GPyOpt, GPy, paramz Building wheel for parameter-sherpa (setup.py) ... done Created wheel for parameter-sherpa: filename=parameter_sherpa-1.0.6-py2.py3-none-any.whl size=542118 sha256=8836fcef51a6e46b135b20539b48b4a7e87ce274b6bb8f630bf39378f80edbfd Stored in directory: /home/fccoelho/.cache/pip/wheels/5e/e2/ac/eb3e175761289cb9a9c5d86045009b9e192ed0df15554f637d Building wheel for GPyOpt (setup.py) ... done Created wheel for GPyOpt: filename=GPyOpt-1.2.6-py3-none-any.whl size=83621 sha256=173cc00fed6f5f8590501dafed53f5f21635365d0961027867aae75d14e8bddb Stored in directory: /home/fccoelho/.cache/pip/wheels/11/b8/44/0282cdff2277bc12f04266de6f104099dec02411879a0ac19f Building wheel for GPy (setup.py) ... done Created wheel for GPy: filename=GPy-1.10.0-cp39-cp39-linux_x86_64.whl size=3286527 sha256=7b5f179a6788cd0a1275f12e328eaaa1309580d1119e93858e863c88422f4bc4 Stored in directory: /home/fccoelho/.cache/pip/wheels/78/fd/57/7c1e4a6f9a5380e2536af9809075ba085b1bb8d38ee84ea183 Building wheel for paramz (setup.py) ... done Created wheel for paramz: filename=paramz-0.9.5-py3-none-any.whl size=102550 sha256=cb0e77e7644d17d8ea46d040c3b03b55de262e17a4fb8106b97b3cabc1efea73 Stored in directory: /home/fccoelho/.cache/pip/wheels/9c/5f/9b/c4273ae8f869387214be2b99598d1b71dbf00672576cb85e74 Successfully built parameter-sherpa GPyOpt GPy paramz Installing collected packages: paramz, MarkupSafe, Werkzeug, Jinja2, itsdangerous, GPy, pymongo, GPyOpt, flask, enum34, parameter-sherpa Attempting uninstall: MarkupSafe Found existing installation: MarkupSafe 1.1.1 Uninstalling MarkupSafe-1.1.1: Successfully uninstalled MarkupSafe-1.1.1 Attempting uninstall: Jinja2 Found existing installation: Jinja2 2.11.2 Uninstalling Jinja2-2.11.2: Successfully uninstalled Jinja2-2.11.2 Successfully installed GPy-1.10.0 GPyOpt-1.2.6 Jinja2-3.0.1 MarkupSafe-2.0.1 Werkzeug-2.0.1 enum34-1.1.10 flask-2.0.1 itsdangerous-2.0.1 parameter-sherpa-1.0.6 paramz-0.9.5 pymongo-3.12.0 WARNING: You are using pip version 21.0.1; however, version 21.2.4 is available. You should consider upgrading via the '/home/fccoelho/Downloads/SageMath/local/bin/python3 -m pip install --upgrade pip' command. 
import sherpa

Começamos por definir o tipo de variável de cada parâmetro, e o tipo de algoritmo de busca que utilizaremos.

parameters = [
    sherpa.Continuous(name='e',range=[0,1]),
    sherpa.Continuous(name='chi',range=[0,.3]),
    sherpa.Continuous(name='phi',range=[0,1]),
    sherpa.Continuous(name='beta',range=[0.1,2]),
    sherpa.Continuous(name='rho',range=[0.2,1]),
    sherpa.Continuous(name='delta',range=[0.01,1]),
    sherpa.Continuous(name='gamma',range=[0.01,1]),
    sherpa.Continuous(name='alpha',range=[0.01,1]),
    sherpa.Continuous(name='mu',range=[0.01,.7]),
    sherpa.Continuous(name='p',range=[0.01,.7]),
    sherpa.Discrete(name='q',range=[1,30]),
    sherpa.Discrete(name='r',range=[1,12]),
    sherpa.Discrete(name='t0',range=[0,25]),
]
algorithm = sherpa.algorithms.RandomSearch(max_num_trials=1000)
# algorithm = sherpa.algorithms.GPyOpt(model_type='GP',max_num_trials=150)

study = sherpa.Study(parameters=parameters,
                     algorithm=algorithm,
                     lower_is_better=True,
                     disable_dashboard=True)

Uma vez definido o problema (study) podemos gerar uma sugestão valores para ver como funciona:

trial = study.get_suggestion()
trial.parameters
{'e': 0.6633981790812898, 'chi': 0.20806102643110655, 'phi': 0.1675576065501, 'beta': 1.6334929266739133, 'rho': 0.4983366877248262, 'delta': 0.3151485461881289, 'gamma': 0.5319575393066565, 'alpha': 0.3685410748547585, 'mu': 0.5881872108774452, 'p': 0.39349886479238544, 'q': 22, 'r': 1, 't0': 19}

Então podemos executar a busca de parâmetros em um simples loop

for trial in study:
    pars = [trial.parameters[n] for n in ['chi', 'phi', 'beta', 'rho', 'delta', 'gamma', 'alpha', 'mu', 'p', 'q', 'r']]
    t0 = trial.parameters['t0']
    T.ode_solve(tspan, inits, num_points=200, params=pars)
    sim = get_sim_array(T.solution)
    H = sim[:ld+t0,-2]*popRJ*trial.parameters['e']
    D = sim[:ld+t0,-1]*popRJ*trial.parameters['e']
    loss = sum((dfRJ.last_available_confirmed-H[t0:t0+ld])**2) +sum((dfRJ.last_available_deaths-D[t0:t0+ld])**2)/2*ld
    study.add_observation(trial=trial,
                          objective=loss,
                          )
    study.finalize(trial)

res = study.get_best_result()
res
{'Trial-ID': 475, 'Iteration': 1, 'alpha': 0.9394810724951352, 'beta': 0.5308498112450182, 'chi': 0.12780938386529786, 'delta': 0.5923550985639527, 'e': 0.9084210097968136, 'gamma': 0.2498206877232359, 'mu': 0.04066936638497329, 'p': 0.3337257898942823, 'phi': 0.09298989032857552, 'q': 21, 'r': 1, 'rho': 0.3594731742120626, 't0': 20, 'Objective': 75443865634.05386}
def plot_results(pars):
    T.ode_solve(tspan, inits, num_points=200, params=list(pars[:-1]))
    t0=pars[-1]
    sim = get_sim_array(T.solution)*popRJ
    h = list_plot(sim[:ld+t0,-2],color='red',legend_label='Cum. cases', plotjoined=True)
    d = list_plot(sim[:ld+t0,-1],color='purple', legend_label='Cum. Deaths', plotjoined=True)
    cc = list_plot(list(zip(range(t0,ld+t0),dfRJ.last_available_confirmed)), color='black',legend_label='cases (obs)')
    cd = list_plot(list(zip(range(t0,ld+t0),dfRJ.last_available_deaths)), color='orange',legend_label='deaths(obs)')
    show(h+d+cc+cd)

plot_results([res['chi'],
              res['phi'],
              res['beta'],
              res['rho'],
              res['delta'],
              res['gamma'],
              res['alpha'],
              res['mu'],
              res['p'],
              res['q'],
              res['r'],
              res['t0']
             ])
Image in a Jupyter notebook

Passo a passo da segunda parte

Resultados

Ter o modelo pronto e adimensionalizado

Encontrar os equilíbrios: ELD (eq. livre de doença) EE (eq. endêmicos)

Por meio de metodo de linearização local, caracterizar as estabilidades dos equilíbrios

Calculo do RO. Analitico e depois substituir os valores para obter uma estimativa numérica.

Gerar simulações e compara com os dados. Descrevam a adequação do modelo para representar a epidemia no país escolhido

Análise de sensibilidade

Otimização dos parâmetros

Estimação bayesiana dos parâmetros.

Discussão

Contextualizar cada resultado

Discutir toda a análise feita de um ponto de vista geral defendendo a aplicabilidade e originalidade dos resultados obtidos.

Implicações do modelo para o controle da COVID

Conclusão

Adicionalmente é necessário escrever um resumo, no início de artigo.

from scipy.interpolate import interp1d
from matplotlib import pyplot as plt

d = [2,6,5,3,8,7,0]
f = interp1d(range(7),d,kind='quadratic')

plt.plot(d, 'o')
t2 = [0, 1, 1.5,2.5,3,3.1,6]
plt.plot(t2,f(t2),'-*' )
[<matplotlib.lines.Line2D object at 0x7f62cf68d2e8>]
Image in a Jupyter notebook
plot(f(t2))
verbose 0 (3797: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 200 points. verbose 0 (3797: plot.py, generate_plot_points) Last error message: ''numpy.ndarray' object is not callable'
/home/fccoelho/Downloads/SageMath/local/lib/python3.7/site-packages/sage/plot/plot.py:2091: DeprecationWarning: elementwise comparison failed; this will raise an error in the future. if funcs == []: /home/fccoelho/Downloads/SageMath/local/lib/python3.7/site-packages/sage/plot/graphics.py:2420: MatplotlibDeprecationWarning: The OldScalarFormatter class was deprecated in Matplotlib 3.3 and will be removed two minor releases later. x_formatter = OldScalarFormatter() /home/fccoelho/Downloads/SageMath/local/lib/python3.7/site-packages/sage/plot/graphics.py:2445: MatplotlibDeprecationWarning: The OldScalarFormatter class was deprecated in Matplotlib 3.3 and will be removed two minor releases later. y_formatter = OldScalarFormatter()
Image in a Jupyter notebook