Welcome to Epidemic Models’ documentation!

Documentation Status GitHub last commit GitHub top language Join discussion Discord

Este projeto tem o objetivo de divulgar como desenvolver modelos para epidemias, desde sua modelagem matemática até sua concepção computacional em Python. Sendo assim composto por vídeos explicativos, notebooks em Python, e diversas visualizações, para ajudar com o entendimento do conteúdo apresentado.


Alguns resultados do século XX…

Para exemplificar o desempenho do conteúdo apresentado, foi feita uma análise com dados do século passado do Reino Unido (United Kingdom), e utilizando os modelos aqui desenvolvidos é possível obter previsões do comportamento dos dados. Essas previsões, juntamente com os dados reais, estão apresentadas na figura a seguir:

Assim como o estudo de correlação entre os modelos SIR obtidos para cada uma das cidades durante os períodos de epidemias do Reino Unido:

_images/UK_models_corr.png

Previsões da COVID-19

Nesta análise utilizamos o modelo desenvolvido para tentar prever os comportamentos da COVID-19. Para isso, primeiramente utilizamos dados de países que já apresentam um comportamento característico da estrutura SIR e já estão no seu período de amortecimento da curva de infectados. Desta forma podemos validar o modelo com relação a sua capacidade de prever eventos futuros, mesmo que somente poucos dias de dados sejam utilizados. Desta forma algumas análises específicas, e de maior impacto, são apresentadas nessa primeira página:

  • Determinação do número básico de reprodução \(R_0\) no decorrer da epidemia
  • Previsão da quantidade de infectados notificados no sistema público
  • Previsão do momento de pico da epidemia

No caso, essas análises foram feitas para Itália 🇮🇹, China 🇨🇳 e Alemanha 🇩🇪, países que já estão em seu período de diminuição do nível de infectados. Para exemplificar a capacidade de previsão do algoritmo de aprendizado desenvolvido, nas visualizações a seguir nós mostramos a previsão feita pelo algoritmo para cada país a medida que os tempo da epidemia foi passando e mais dados foram utilizados para a aprendizagem:

Nosso modelo aprende com algoritmos de otimização os parâmetros da estrutura SIR, juntamente com a proporção da população que está sendo registrada pelo sistema de saúde. No caso, ele aprende três parâmetros fundamentais do modelo: \(\beta\) (contatos por dia), \(r\) (em que \(1/r\) é o tempo médio de recuperação da doença), e o \(S(0)\) (quantidade inicial de suscetíveis), para que a epidemia tenha o comportamento que os dados mostram. No caso os valores ajustados para cada país foram:

País

\(\beta\)

\(r\)

Itália 🇮🇹

0.2038

0.0233

China 🇨🇳

0.3133

0.0445

Alemanha 🇩🇪

0.2195

0.0576

Caso queira tentar ajustar os parâmetros você mesmo, basta clicar no botão abaixo!


Determinação do \(R_0\)

Note que o parâmetro \(R_0\), é determinado a partir dos dois outros parâmetros característicos do modelo SIR, \(R_0 = \beta / r\). Aqui utilizamos um modelo SIR que pondera a quantidade de suscetíveis, uma vez que nem toda a população pode ser considerada suscetível, visto que nem todas as pessoas infectadas e recuperadas são notificadas ao sistema público. Além disso, nem toda a população é exposta ao vírus, devido a políticas públicas, isolamentos, etc.

Note

Desta forma, note que o modelo desenvolvido somente modela as pessoas notificadas pelo sistema de saúde, sendo assim, representativo de uma parte da verdade situação do país.

Dito isso, é possível definir os valores encontrados pelo algoritmo de aprendizagem para o parâmetro \(R_0\) a medida que os dias da pandemia passaram:

Previsões do consumo do sistema público

Um dos parâmetros que nosso algoritmo aprende durante seu processo de treinamento é um parâmetro que pondera a quantidade da população de suscetíveis inicial (simplesmente uma técnica para melhorar o condicionamento númerico do algoritmo). Porém, como esse parâmetro tende sempre a estimar o valor de \(S(0)\) igual ao valor de \(R(\infty)\). Note que sempre é verdade \(S(0) \geq R(\infty)\). Como os dados medidos são somente das pessoas notificadas e acompanhadas pelo sistema de saúde, podemos concluir que o valor de \(R(\infty)\) é a quantidade de pessoas que foram contaminadas, frequentaram o sistema de saúde para o diagnóstico e por isso estão na base de dados. Nosso algoritmo prevê a quantidade de \(S(0) = R(\infty)\), desta forma para cada novo dia de dados temos uma nova previsão de qual será a quantidade de pessoas absorvidas pelo sistema de saúde. Nos gráficos a seguir conseguimos mostrar o erro percentual entre o real valor de \(R(\infty)\) e o valor estimado por nosso modelo a cada dia da epidemia:

Previsões dos picos epidêmicos da COVID-19

Nesta análise apresentamos o efeito da quantidade de dados na performance do modelo desenvolvido analisando a capacidade de prever o dia em que acontecerá o pico da quantidade de infectados da epidemia da COVID-19. Para isso estamos utilizando dados de países que já tiveram seu pico de contágio e atualmente estão no período de amortecimento da quantidade de infectados. Alguns dos países analisados foram a China 🇨🇳, Itália 🇮🇹 e Alemanha 🇩🇪, que possibilitaram as análises abaixo. Nestas figuras é mostrado o erro do modelo ao tentar prever o dia de pico, para cada dia da epidemia:

Autores

Modelo SIR

[1]:

from IPython.display import IFrame

# Youtube
IFrame("https://www.youtube.com/embed/v-91XHpjugo", "100%", 500)

[1]:

Configurações

[2]:
import numpy as np
import pandas as pd
from scipy.integrate import odeint

from bokeh.palettes import brewer
from bokeh.plotting import figure, show
from bokeh.io import output_notebook

output_notebook()
Loading BokehJS ...

Suscetível, Infectado e Removido (Recuperado)

A proposta é introduzir um modelo determinístico para análise de evolução de uma epidemia. O modelo aqui representado foi inicialmente proposto por Kermack, W. and McKendrick, A., 1927. A contribution to the mathematical theory of epidemics. Proc. R. Soc. London A 115, 700-721. O modelo tem as seguintes premissas:

  1. Uma posulação invariante no tempo com N individuos (fechada).
  2. Taxas constantes: contato e remoção.
  3. Desconsiderando variações demográficas: nascimentos e mortes.
  4. População “bem misturada”.

O modelo é fundamentado em três estados: S - suscetível, I - infectado e R - removido ou recuperado.

\[S(t)+I(t)+R(t)= N\]
[3]:

# Tamanho da populção - N
N = 500

# Valores iniciais
I0 = 1
R0 = 0
S0 = N - I0

O conjunto de equações diferenciais que caracteriza o modelo é descrito abaixo. No modelo

\[\beta - \text{representa a taxa de transmissão ou taxa efetiva de contato}\]

e

\[r - \text{a taxa de remoção ou recuperação.}\]
\[\begin{split}\begin{split} \frac{dS(t)}{dt} & = -\beta S(t) I(t) \\ \frac{dI(t)}{dt} & = \beta S(t) I(t) - rI(t) \\ \frac{dR(t)}{dt} & = r I(t) \end{split}\end{split}\]
[4]:

# Definição do conjunto de equações diferencias não lineares que formam o modelo.

def SIRm(y, t, N, Beta, r):
    S, I, R = y
    Sdot = -(Beta * S * I)
    Idot = (Beta * S * I)  - r * I
    Rdot = r * I
    return Sdot, Idot, Rdot

Analisando a evolução de \(\beta\) e \(r\)

[5]:

# Beta -  taxa de contato,
# r - taxa média de recuperação (in 1/dia).
Beta = 0.0009
r = 1e-1

# Resolução da simulação - Escala temporal (dias)
t = np.linspace(0, 100, 1000)

# Vetor de condições iniciais
y0 = S0, I0, R0

# Integrando as Equações do modelo SIR
ret = odeint(SIRm, y0, t, args=(N, Beta, r))
S, I, R = ret.T

# Visualizando a evolução da Epidemia - S(t), I(t) e R(t)
p = figure(tools="pan,hover,lasso_select",
           x_range=(0, t[-1]), y_range=(0, N),
           plot_width=600, plot_height=400)

p.line(t, S, legend_label="Suscetíveis", color="#ffd885", line_width=3)
p.line(t, I, legend_label="Infectados",  color="#de425b", line_width=3)
p.line(t, R, legend_label="Removidos",   color="#99d594", line_width=3)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"
p.legend.click_policy = "mute"
p.legend.items.reverse()

show(p)

Salvando o modelo

[6]:

import pandas as pd

output = pd.DataFrame({"S": S, "I": I, "R": R})
output.to_csv('SIRpd', index=False)

print(output)

              S         I           R
0    499.000000  1.000000    0.000000
1    498.954252  1.035561    0.010187
2    498.906883  1.072381    0.020736
3    498.857833  1.110506    0.031660
4    498.807046  1.149981    0.042973
..          ...       ...         ...
995    5.855617  0.232760  493.911622
996    5.855495  0.230564  493.913941
997    5.855374  0.228388  493.916238
998    5.855254  0.226232  493.918514
999    5.855135  0.224097  493.920768

[1000 rows x 3 columns]

Modelo Reed Frost

[1]:

from IPython.display import IFrame

# Youtube
IFrame("https://www.youtube.com/embed/D8S4GGmHUMw", "100%", 500)

[1]:

Configurações

[2]:

import numpy as np
from scipy.integrate import odeint

from bokeh.palettes import brewer, Inferno10
from bokeh.plotting import figure, show
from bokeh.io import output_notebook
from bokeh.layouts import gridplot

output_notebook()

Loading BokehJS ...

Modelo para pequeno grupo

Há no senso comum o entendimento que há uma chance de “pegar” uma doença. Supondo que a probailidade que representa esta sensação seja \(\mathbf{p}\) e que a probilidade de escapar desta doença seja \(\mathbf{q}\). Vamos admitir um cenário onde tenhamos uma família de três pessoas \(\mathbf{N}=3\) e um dos individuos estaje infectado (I). Supondo que neste caso \(\mathbf{p}=0.4\) logo \(\mathbf{q}=0.6\). Admitindo este cenário pode-se explorar as seguintes progressões da doença (\(P_i\) - infectado e \(\overline{P_i}\) - não infectado) :

Pensando em casos:

  1. Apenas o primeiro ifectado desenvovle a doença.
  2. O infectado transmite para um dos outros dois suscetíveis apenas.
  3. Os dois sucetíveis são infectados um após o outro.
  4. Ambos os suscetíveis são infectados simultaneamente.

Traduzindo em cadeias de transmissão:

  1. 1 \(\rightarrow\) 0
  2. 1 \(\rightarrow\) 1 \(\rightarrow\) 0
  3. 1 \(\rightarrow\) 1 \(\rightarrow\) 1
  4. 1 \(\rightarrow\) 2

Traduzindo em probabilidades de cada cadeia de transmissão:

  1. \(\overline{P_1}\) e \(\overline{P_2}\) - \(\mathbf{q^2}\)
  2. \(\overline{P_1}\) e \(P_2\) - \(\mathbf{pq^2}\) ou \(P_1\) e \(\overline{P_2}\) - \(\mathbf{pq^2}\)
  3. \(\overline{P_1}\) e \(P_2\) \(\rightarrow\) \(P_1\) e \(P_2\) - \(\mathbf{p^2q}\) ou \(P_1\) e \(\overline{P_2}\) \(\rightarrow\) \(P_1\) e \(P_2\) - \(\mathbf{p^2q}\)
  4. \(P_1\) e \(P_2\) - \(\mathbf{p^2}\)
[3]:

p = 0.4
q = 1 - p

C1 = q**2
C2 = 2*p*(q**2)
C3 = 2*q*(p**2)
C4 = p**2

print(C1,C2,C3,C4, C3+C4, C1+C2+C3+C4)

0.36 0.288 0.19200000000000003 0.16000000000000003 0.3520000000000001 1.0

Considerando as probabilidades calculadas para cada uma das sequências, aplicada um cenário de 1000 famílias com três integrantes, pode-se estimar que cerca de 360 famílias, atigidas pela doença, terão apenas 1 infectado, 288 famílias terão pelo menos 1 infectado e 352 famílias todos serão infectados.

Se generalizarmos a análise pode-se dizer que o número de infectados \(\mathbf{I_{t}}\) na próxima geração de infectados será \(\mathbf{I_{t+1}}\). O termo \(\mathbf{t}\) indica a geração da epidemia e indica com uma certa imprecisão a escala de tempo.

\begin{align*} P(I_{t+1}=i_{t+1}\vert S_t=s_t,I_t=i_t) \sim \binom{s_t}{i_t} (1-q^i_t)^{I_{t+1}} (q)^{i_t(s_t-i_{t+1})}, s_t \geq i_{t+1} \end{align*}

Generalizando o modelo

A partir da proposta que o contágio é dependente de uma probabilidade \(\mathbf{p}\) de transmissão da doença, que corresponde a uma probabilidade \(\mathbf{q}\) de escapar de ser contaminado, pode-se estabelecer um modelo que seja capaz de representar a cadeia de transmissão.

\[\begin{split}\begin{split} I_{t+1} & \sim binomial(S(t),p) \\ S_{t+1} & = S_t - I_{t+1} \\ R_{t+1} & = R_t + I_t \end{split}\end{split}\]

Observe a semelhança com o modelo determinístico, exceto o fato do número de infectados ser um valor randômico, governado pela distribuição Binomial.

[4]:

# Modelo Reed-Frost
ngen  = 30;     # número de gerações
Sinit = 2000;   # população suscetível
Iinit = 1;      # infectados
Rinit = 0;      # removidos / recuperados

q = 0.999;

nsims = 2000;   # número de simulações

x, y = list(), list()
M = np.zeros((nsims-1,ngen-1))

TOOLS="zoom_in,zoom_out,save"
p = figure(tools=TOOLS, plot_width=600, plot_height=400)

for i in range(1,nsims-1):
    S = np.linspace(0,1,ngen) - np.linspace(0,1,ngen)
    I = np.linspace(0,1,ngen) - np.linspace(0,1,ngen)
    R = np.linspace(0,1,ngen) - np.linspace(0,1,ngen)
    S[1], I[1], R[1] = Sinit, Iinit, Rinit
    for j in range(1,ngen-1):
        # np.random.binomial(n, p, 1000)
        I[j+1] = np.random.binomial(S[j], 1 - q**I[j], 1);
        S[j+1] = S[j] - I[j+1]
        R[j+1] = R[j] + I[j]
        M[i,j] = I[j]

    p.line(range(ngen), I, line_width=2, color="#8e44ad", line_alpha=0.05)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Infectados"
p.xaxis.axis_label = "Gerações"

show(p)

Análise dos cenários para cada geração nova de infectados

Observe que para cada nova geração a distribuição estatística varia e pode-se observar o comportamento da distribuição Binomial. Neste sentido pode-se observar quem “governa” a incerteza da evolução da epidemia. A linha pontilhada em cada histograma representa a média para cada geração.

[5]:

num = 10
x = np.linspace(0, 100, 1000)
pallete = Inferno10

fig_list = []
for i in range(num):
    hist, edges = np.histogram(M[:,i], density=False, bins=50)
    p = figure(tools="pan,hover,lasso_select", plot_width=300, plot_height=300)
    p.quad(top=hist, bottom=0, left=edges[:-1], right=edges[1:],
     fill_color=pallete[i], fill_alpha=0.5, line_color=pallete[i],
     legend_label="k = " + str(i) )
    p.grid.grid_line_alpha = 0
    p.ygrid.band_fill_alpha = 0.1
    p.ygrid.band_fill_color = "olive"
    p.yaxis.axis_label = "Frequência"
    p.xaxis.axis_label = "Infectados"
    fig_list.append(p)

grid = gridplot(fig_list, ncols=2)
show(grid)

[6]:

num = 15

p = figure(tools="pan,hover,lasso_select", y_axis_type="log", plot_width=600, plot_height=400)

for i in range(num):
    x = (i+1) * np.ones((len(M[:,i]),))
    try:
        p.scatter(x, M[:,i], size=5, fill_color=pallete[i], fill_alpha=0.05, line_alpha=0)
    except:
        p.scatter(x, M[:,i], size=5, fill_color=pallete[i-10], fill_alpha=0.05, line_alpha=0)
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_alpha = 0.1
p.ygrid.band_fill_color = "olive"
p.yaxis.axis_label = "Infectados"
p.xaxis.axis_label = "Gerações"

show(p)

Evolução da média e desvio padrão

Para uma visão mais clara podemos visualizar a média e o desvio padrão para cada geração.

[7]:
a = M.mean(0)
e = M.std(0)
x = np.linspace(0,ngen, ngen-1)

p = figure(tools="pan,hover,lasso_select", plot_width=600, plot_height=400)

p.varea(x=x, y1=a+e/2, y2=a-e/2, legend_label="Desvio Padrão", fill_alpha=0.3, fill_color="#3498db")
p.line(x, a, legend_label='Média', line_cap="round", color='#e67e22', line_width=4)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_alpha = 0.1
p.ygrid.band_fill_color = "olive"
p.yaxis.axis_label = "Infectados"
p.xaxis.axis_label = "Geração"

show(p)

SIR Dados Artificiais

[1]:
import pylab as pp
import numpy as np
import pandas as pd

from scipy import integrate, interpolate, optimize
from scipy.integrate import odeint

from bokeh.palettes import brewer
from bokeh.plotting import figure, show
from bokeh.io       import output_notebook

output_notebook()
Loading BokehJS ...

Dados artificiais

Uma prática interessante na análise de dados é testarmos se o algoritmo proposto e principalmente sua estrutura está consistente. Podemos então utilizar dados artificiais para avaliar se o processo de ajuste está estruturado corretamente.

[2]:

# Lendo o arquivo de dados no formato 'filename.csv'
data = pd.read_csv("./PG_IMT/DadosEpidemia/SIRpd")
# Preview das cinco primeiras linhas
data.head()

[2]:
S I R
0 499.000000 1.000000 0.000000
1 498.949042 1.040746 0.010213
2 498.896013 1.083146 0.020841
3 498.840830 1.127267 0.031903
4 498.783405 1.173179 0.043415
[3]:

s_array = data[["S", "I", "R"]].to_numpy()

Sd = s_array[:,0]
Id = s_array[:,1]
Rd = s_array[:,2]

Gerando ruído gaussiano

[4]:

Sdn = np.random.normal(0, np.mean(Sd)/10, len(Sd))
Idn = np.random.normal(0, np.mean(Id)/5,  len(Sd))
Rdn = np.random.normal(0, np.mean(Rd)/10, len(Sd))

Sd = Sd + Sdn
Id = Id + Idn
Rd = Rd + Rdn

Visualizando os dados

[5]:

# Visualizando a evolução da Epidemia - S(t), I(t) e R(t)
t = np.linspace(0,len(Sd),len(Sd))

TOOLS="hover,crosshair,pan,wheel_zoom,zoom_in,zoom_out,box_zoom,undo,redo,reset,tap,save,box_select,poly_select,lasso_select,"

p = figure(tools=TOOLS, plot_width=600, plot_height=400)

p.scatter(t, Sd, legend_label="Suscetíveis - dados",
          radius=3.8, fill_color="#ffd885", fill_alpha=0.6, muted_color="#ffd885", muted_alpha=0.2, line_color=None)
p.scatter(t, Id, legend_label="Infectados - dados",
          radius=3.8, fill_color="#de425b", fill_alpha=0.6, muted_color="#de425b", muted_alpha=0.2, line_color=None)
p.scatter(t, Rd, legend_label="Removidos - dados",
          radius=3.8, fill_color="#99d594", fill_alpha=0.6, muted_color="#99d594", muted_alpha=0.2, line_color=None)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"
p.legend.click_policy = "mute"
p.legend.items.reverse()

show(p)

O problema

O conjunto de equações diferenciais que caracteriza o modelo é descrito abaixo. No modelo $:nbsphinx-math:beta `- :nbsphinx-math:text{representa a taxa de transmissão ou taxa efetiva de contato}` $ e \(r - \text{a taxa de remoção ou recuperação.}\)

\[\begin{split}\begin{split} \frac{dS(t)}{dt} & = -\beta S(t) I(t) \\ \frac{dI(t)}{dt} & = \beta S(t) I(t) - rI(t) \\ \frac{dR(t)}{dt} & = r I(t) \end{split}\end{split}\]

Gostaríamos de identificar quais parâmetros \(\beta\) e \(r\) resultam num melhor ajuste do modelo para os dados de S,I e R

[6]:

def SIRmodel(y, t, Beta,r):
    S, I, R = y
    Sdot = -(Beta * S * I)
    Idot = (Beta * S * I)  - r * I
    Rdot = r * I
    return Sdot, Idot, Rdot

# Resolução da simulação - Escala temporal (dias)

Obtendo \(y_s(\theta,k) = [S \; I \: R]\)

O trecho a seguir retorna os valores sintetizados \(y_s(\theta,k) = [S \; I \: R]\) representa o dado sintetizado a partir de um modelo sintetizado para uma determinada amostra \(k\) e \(\theta\) representa o vetor ed parâmetros \(\theta = [ \beta \; \; r]^T\). A partir de uma condição inicial \(y_0\).

[7]:

def SIRsim(y0,t,theta):
    Beta = theta[0]
    r = theta[1]
    ret = integrate.odeint(SIRmodel,y0,t,args=(Beta,r))
    S, I, R = ret.T
    return S, I, R

Condições inicias - \(y_0\) e \(\theta_0\)

[8]:

# Tamanho da populção - N
N = 500

# Valores iniciais
I0, R0 = 1, 0
S0 = N - I0

# Vetor de condições iniciais

y0 = S0, I0, R0

# Beta -  taxa de contato,
# r - taxa média de recuperação (in 1/dia).

theta0 = [1e-4,1e-2] # valores iniciais

# Definição do conjunto de equações diferencias não lineares que formam o modelo.

t = np.linspace(0, 1000, 1000)

Estimação de parâmetros

Para estimarmos os parâmetros do modelo \(\mathbf{\beta}\) e \(\mathbf{r}\), vamos utilizar inicialmente o método de mínimos quadrados. Podemos então formular o problema a partir da Equação abaixo. Na Equação \(y_m(k)\) representa o dado real em cada amostra \(k\); \(y_s(\theta,k)\) representa o valor estimado a partir da simulação do modelo para uma determinada amostra \(k\) e \(\theta\) representa o vetor ed parâmetros \(\theta = [ \beta \; \; r]^T\).

\[min_{\theta}= \sum_{k=1}^{K}(y_m(k) - y_s(\theta,k))^2\]

A equação formula a pergunta: quais os valores de \(beta\) e \(r\) que minizam o erro quadrático quando comparados com os dados reais.

[9]:

def ErroQuadratico(Sd,Id,Rd,y0,t,theta0):
    """ function to pass to optimize.leastsq
        The routine will square and sum the values returned by
        this function"""
    [S,I,R] = SIRsim(y0,t,theta0)
    erroS = S - Sd
    erroI = I - Id
    erroR = R - Rd
    EQ = np.concatenate([erroI,erroR])
    return EQ

def objetivo(p):
    return ErroQuadratico(Sd,Id,Rd,y0,t,p)

Minimização da função custo

[10]:

(c,kvg) = optimize.leastsq(objetivo,theta0)
print(c)

[9.93827391e-05 9.96628326e-03]

Visualização

[11]:

[Sa,Ia,Ra] = SIRsim(y0,t,c)

# Visualizando a evolução da Epidemia - S(t), I(t) e R(t)
p.line(t, Sa, legend_label="Suscetíveis", color="#f9bd3d", line_width=3)
p.line(t, Ia, legend_label="Infectados",  color="#f4193c", line_width=3)
p.line(t, Ra, legend_label="Removidos",   color="#45c83a", line_width=3)

show(p)

SIR 🇬🇧 Reino Unido 1944-1964

[1]:
import pylab as pp
import numpy as np
import pandas as pd
import scipy.signal as scs

from scipy import integrate, interpolate, optimize
from scipy.integrate import odeint
from pylab import *
import netCDF4

from bokeh.models   import ColumnDataSource, RangeTool, LinearAxis, Range1d
from bokeh.palettes import brewer, Inferno10
from bokeh.plotting import figure, show
from bokeh.layouts  import column
from bokeh.io       import output_notebook

output_notebook()
Loading BokehJS ...

Dados Reais

Dados referentes a cidade de Londres entre 1944 e 1964. População de aproximadamente 2.5 milhões de habitantes na época.

[2]:

# Lendo o arquivo de dados no formato 'filename.csv'
data = pd.read_csv("./PG_IMT/DadosEpidemia/UKCities/London.csv")
# Preview das cinco primeiras linhas
data.head()

s_array = data[["cases", "births", "pop","time"]].to_numpy()

Id = s_array[:,0]
Bd = s_array[:,1]
Sd = s_array[:,2]
t  = s_array[:,3]

Visualizando os dados

[3]:

TOOLS="zoom_in,zoom_out,save"
p1 = figure(tools=TOOLS, plot_width=600, plot_height=300)
p2 = figure(tools=TOOLS, plot_width=600, plot_height=300)
p3 = figure(tools=TOOLS, plot_width=600, plot_height=300)

p1.line(data["time"], data["births"], line_width=4, color="#8e44ad", line_alpha=0.9)

p1.grid.grid_line_alpha = 0
p1.ygrid.band_fill_color = "olive"
p1.ygrid.band_fill_alpha = 0.1
p1.yaxis.axis_label = "Nascimentos"
p1.xaxis.axis_label = "Data"

p2.line(data["time"], data["pop"], line_width=4, color="#8e44ad", line_alpha=0.9)

p2.grid.grid_line_alpha = 0
p2.ygrid.band_fill_color = "olive"
p2.ygrid.band_fill_alpha = 0.1
p2.yaxis.axis_label = "População"
p2.xaxis.axis_label = "Data"

p3.line(data["time"], data["cases"], line_width=4, color="#8e44ad", line_alpha=0.9)

p3.grid.grid_line_alpha = 0
p3.ygrid.band_fill_color = "olive"
p3.ygrid.band_fill_alpha = 0.1
p3.yaxis.axis_label = "Casos"
p3.xaxis.axis_label = "Data"

show(column(p1,p2,p3))

Selecionando um ano específico para análise

[4]:

selected_year = 1948

is_1948 = (data['time'].astype(int) == selected_year)
LD48 = data[is_1948]
LD48.head()

p = figure(tools=TOOLS, plot_width=600, plot_height=300)
p.line(LD48["time"], LD48["cases"], line_width=4, color="#8e44ad", line_alpha=0.9)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Casos"
p.xaxis.axis_label = "Data"

show(p)

Os dados

Conhecemos o número de casos, os nascimentos e a população total. Podemos determinar os demais valores a partir da relação \(S(t)+I(t)+R(t)=N\) onde N vamos considerar a população.

[5]:

s_array = LD48[["cases", "births", "pop","time"]].to_numpy()

Id = s_array[:,0]
Bd = s_array[:,1]
Sd = s_array[:,2] - Id + Bd
t  = s_array[:,3]

O problema

O conjunto de equações diferenciais que caracteriza o modelo é descrito abaixo. No modelo \(\beta\) - representa a taxa de transmissão ou taxa efetiva de contato e \(r\) - a taxa de remoção ou recuperação.

\[\begin{split}\begin{align} \frac{dS(t)}{dt} & = -\beta S(t) I(t) \\ \frac{dI(t)}{dt} & = \beta S(t) I(t) - rI(t) \\ \frac{dR(t)}{dt} & = r I(t) \end{align}\end{split}\]

Gostaríamos de identificar quais parâmetros \(\beta\) e \(r\) resultam num melhor ajuste do modelo para os dados de S, I e R

[6]:

def SIRmodel(y, t, Beta, r):
    S, I = y
    Sdot = -Beta * S * I
    Idot = Beta * S * I  - r * I
    return Sdot, Idot

# Resolução da simulação - Escala temporal (dias)

Obtendo \(y_s(\theta,k) = [S \; I \: R]\)

O trecho a seguir retorna os valores sintetizados \(y_s(\theta,k) = [S \; I \: R]\) representa o dado sintetizado a partir de um modelo sintetizado para uma determinada amostra \(k\) e \(\theta\) representa o vetor ed parâmetros \(\theta = [ \beta \; \; r]^T\). A partir de uma condição inicial \(y_0\).

[7]:

def SIRsim(y0, t, theta):
    Beta, r = theta[0], theta[1]
    ret = integrate.odeint(SIRmodel, y0, t, args=(Beta,r))
    S, I = ret.T
    return S, I

Condições inicias - \(y_0\) e \(\theta_0\)

[8]:

# Valores iniciais
I0 = Id[0]
S0 = Sd[0]

# Vetor de condições iniciais

y0 = S0, I0

# Beta -  taxa de contato
# r - taxa média de recuperação (in 1/dia).

# theta = [Beta, r]
theta0 = [1e-8, 1e-1] # valores iniciais

# Definição do conjunto de equações diferencias não lineares que formam o modelo.
t = (t - 1948) * 365

Estimação de parâmetros

Para estimarmos os parâmetros do modelo \(\mathbf{\beta}\) e \(\mathbf{r}\), vamos utilizar inicialmente o método de mínimos quadrados. Podemos então formular o problema a partir da Equação abaixo. Na Equação \(y_m(k)\) representa o dado real em cada amostra \(k\); \(y_s(\theta,k)\) representa o valor estimado a partir da simulação do modelo para uma determinada amostra \(k\) e \(\theta\) representa o vetor ed parâmetros \(\theta = [ \beta \; \; r]^T\).

\[min_{\theta}= \sum_{k=1}^{K}\left(y_m(k) - y_s(\theta,k) \right)^2\]

A equação formula a pergunta: quais os valores de \(beta\) e \(r\) que minizam o erro quadrático quando comparados com os dados reais.

[9]:

def ErroQuadratico(Sd,Id,y0,t,theta0,w):
    """ function to pass to optimize.leastsq
        The routine will square and sum the values returned by
        this function"""
    [S,I] = SIRsim(y0,t,theta0)
    erroS = w[0] * (S - Sd)
    erroI = w[1] * (I - Id)
    EQ = np.concatenate([erroS,erroI])
    if sum(np.isnan(EQ)) >= 1:
        print("Error!!!")
    return EQ

def objetivo(p, S, I, y0, t, w):
    return ErroQuadratico(S,I,y0,t,p,w)

Minimização da função custo

[10]:

# Criação das ponderações dos erros
wS = max(Id) / max(Sd)
w = [wS, 1]

(c,kvg) = optimize.leastsq(objetivo, theta0, args=(Sd, Id, y0, t, w))
print(c)

[2.10765916e-08 6.83427510e-02]

Visualização

[11]:

[Sa,Ia] = SIRsim(y0,t,c)

p = figure(tools=TOOLS, y_range=(0,5000), plot_width=600, plot_height=400)

p.scatter(t, Sd, legend_label="Suscetíveis - dados", size=8, fill_color="#ffd885", fill_alpha=0.7, line_alpha=0)
p.scatter(t, Id, legend_label="Infectados - dados", size=8, fill_color="#de425b", fill_alpha=0.7, line_alpha=0)

p.line(t, Sa, legend_label="Suscetíveis", line_width=4, color="#f9bd3d", line_cap='round', line_alpha=0.9)
p.line(t, Ia, legend_label="Infectados", line_width=4, color="#f4193c", line_cap='round', line_alpha=0.9)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

show(p)

Reamostrando os dados

[12]:

Sd_res, t_res = scs.resample(Sd, 365, t=t)
Id_res, t_res = scs.resample(Id, 365, t=t)

p = figure(tools=TOOLS, y_range=(0,4000), plot_width=600, plot_height=400)

p.scatter(t, Id, legend_label="Infectados", size=10, fill_color="#f4193c", fill_alpha=0.6, line_alpha=0)
p.scatter(t_res, Id_res, legend_label="Infectados - Resample", size=4, fill_color="#8e44ad", line_alpha=0)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

show(p)

Restimando os parâmetros

[13]:

(c, kvg) = optimize.leastsq(objetivo, theta0, args=(Sd_res, Id_res, y0, t_res, w))
print(c)

[Sa,Ia] = SIRsim(y0, t_res, c)

p.line(t_res, Ia, legend_label="Infectados - Modelo", line_width=4, color="#f4193c", line_cap='round', line_alpha=0.9)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

show(p)

[2.15762357e-07 6.89608748e-01]

Outros métodos de estimação

[14]:

import warnings
warnings.filterwarnings("error")

from scipy.optimize import dual_annealing, shgo, basinhopping, differential_evolution

def cost_function(theta0,Sd,Id,y0,t,w):
    """ function to pass to optimize.leastsq
        The routine will square and sum the values returned by
        this function"""
    try:
        [S, I] = SIRsim(y0, t, theta0)
        erroS = (w[0]*(S - Sd)**2)
        erroI = (w[1]*(I - Id)**2)
        final_error = sum(erroI)
    except:
        final_error = 10**14
    return final_error


beta_approx = 1 / data["pop"].max()
r_approx = 1 / 365

x0 = [beta_approx, r_approx]
lower, upper = [x0[0]/10, x0[1]/10], [10*x0[0], 1000*x0[1]]
bounds = zip(lower, upper)

#res = dual_annealing(cost_function, x0=x0, maxiter=2000, bounds=list(bounds), args=(Sd_res, Id_res, y0, t_res, w))
#res = shgo(cost_function, bounds=list(bounds), args=(Sd_res, Id_res, y0, t_res, w))
#res = basinhopping(cost_function, x0=x0, niter=200, minimizer_kwargs=dict(args=(Sd_res, Id_res, y0, t_res, w)))

res = differential_evolution(cost_function, list(bounds), args=(Sd_res, Id_res, y0, t_res, w))
res

[14]:
     fun: 43830917.444261
     jac: array([2.89188102e+17, 3.06168899e+07])
 message: 'Optimization terminated successfully.'
    nfev: 906
     nit: 26
 success: True
       x: array([2.16790744e-07, 6.92576449e-01])
[15]:

[Sa,Ia] = SIRsim(y0, t_res, res.x)

p = figure(tools=TOOLS, y_range=(0,4000), plot_width=600, plot_height=400)

p.scatter(t, Id, legend_label="Infectados", size=10, fill_color="#f4193c", fill_alpha=0.6, line_alpha=0)
p.scatter(t_res, Id_res, legend_label="Infectados - Resample", size=4, fill_color="#8e44ad", line_alpha=0)

p.line(t_res, Ia, legend_label="Infectados - Modelo", line_width=4, color="#f4193c", line_cap='round', line_alpha=0.9)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

show(p)

Análise anual

[16]:

years = data["time"].astype(int).unique().tolist()

p = figure(tools=TOOLS+",hover",
            x_range=(min(years), max(years)),
            plot_width=600, plot_height=300)

p.line(data["time"], data["cases"],
       legend_label="Casos", line_width=4, color="#f4511e", line_cap='round', line_alpha=0.9)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"
p.toolbar.autohide = True

show(p)

Detecção dos períodos de contágio

[17]:

from PyAstronomy import pyasl

# Filtrando os dados para o processo de derivação
filt_cases = pyasl.smooth(data["cases"], 13, 'hamming')


# Incluindo os dados filtrados no plot
p.line(data["time"], filt_cases, line_dash="dotted",
       legend_label="Casos Filtrado", line_width=3, color="#8e44ad", line_cap='round', line_alpha=0.9)

show(p)

[18]:

cases_variation = np.diff(filt_cases) # Calculando a taxa de variação

threshold = np.std(cases_variation)   # Calculando o threshold baseado no desvio
                                      # padrão da taxa de variação
threshold_vec = [threshold for k in cases_variation]

p = figure(tools=TOOLS+",hover",
            plot_width=600, plot_height=300)

p.line(data["time"][:-1], cases_variation,
       legend_label="Derivada dos casos", line_width=4, color="#f4511e", line_cap='round', line_alpha=0.9)
p.line(data["time"][:-1], threshold_vec,
       legend_label="Threshold", line_width=4, color="#8e44ad", line_cap='round', line_alpha=0.9)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"
p.toolbar.autohide = True

show(p)

[19]:

def findEpidemyBreaks(cases, threshold_prop, cases_before=10):
    """
    """

    # Filtering the data
    filt_cases = pyasl.smooth(cases, 11, 'hamming')
    # Compute the derivative and standard deviation
    cases_variation = np.diff(filt_cases).tolist()
    threshold = threshold_prop * np.std(cases_variation)

    start_points, end_points = [], []
    in_epidemy, out_epidemy = False, False
    for k, value in enumerate(cases_variation):
        if not in_epidemy:
            # Check value
            if value > threshold:
                in_epidemy = True
                # Find the start point
                start_index = 0 if k-cases_before < 0 else k-cases_before
                window = [abs(v) for v in cases_variation[start_index:k]]
                ref_index = window.index(min(window))
                start_points.append(k - (cases_before - ref_index))
        else:
            check_1 = (cases_variation[k-1] < 0)
            check_2 = (value >= 0)
            if check_1 and check_2:
                in_epidemy = False
                end_points.append(k)
    return start_points, end_points

[20]:

start, end = findEpidemyBreaks(data["cases"], 1, 12)

p = figure(tools=TOOLS+",hover",
            x_range=(min(years), max(years)),
            y_range=(0, 8000),
            plot_width=600, plot_height=400)

p.line(data["time"], data["cases"],
       legend_label="Casos", line_width=4, color="#f4511e", line_cap='round', line_alpha=0.9)



windowed_data = { "I":[], "S":[], "B":[], "t":[] }
for s, e in zip(start, end):
    windowed_data["B"].append(data["births"][s:e].to_numpy())
    windowed_data["S"].append(data["pop"][s:e].to_numpy())
    windowed_data["I"].append(data["cases"][s:e].to_numpy())
    windowed_data["t"].append(data["time"][s:e].to_numpy())

    p.line(windowed_data["t"][-1], windowed_data["I"][-1],
          legend_label="Casos na Epidemia", line_width=3, color="#8e44ad", line_cap='round', line_alpha=0.9)

    p.line([data["time"][s], data["time"][s]], [0, 10000], line_dash="dotted",
          legend_label="Início da Epidemia", line_width=1, color="#455a64", line_cap='round', line_alpha=0.9)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Ano"
p.toolbar.autohide = True


index = 5
days = (windowed_data["t"][index] - windowed_data["t"][index][0]) * 365

p1 = figure(tools=TOOLS+",hover",
            plot_width=600, plot_height=400)

p1.line(days, windowed_data["I"][index],
       legend_label="Casos", line_width=4, color="#8e44ad", line_cap='round', line_alpha=0.9)

p1.grid.grid_line_alpha = 0
p1.ygrid.band_fill_color = "olive"
p1.ygrid.band_fill_alpha = 0.1
p1.yaxis.axis_label = "Indivíduos"
p1.xaxis.axis_label = "Dias"
p1.toolbar.autohide = True

show(column(p,p1))

[21]:

final_data = {
    "data": {
        "original": [],
        "resampled": [],
        "simulated": []
    },
    "pars": {
        "beta": [],
        "r": []
    }
}

zipped_data = zip(
    windowed_data["S"],
    windowed_data["I"],
    windowed_data["B"],
    windowed_data["t"])

for S, I, B, t in zipped_data:
    year_ref = t[0]
    S = S - I + B # Calculando os suscetiveis
    t = (t - year_ref) * 365 # Tempo em dias
    # Condições iniciais
    y0 = S[0], I[0]
    theta0 = [1e-8, 1e-1]
    # Criando a ponderação
    # dos erros
    wS = max(I) / max(S)
    w = [wS, 1]
    # Reamostrando os dados
    Sd_res, t_res = scs.resample(S, int(t[-1]), t=t)
    Id_res, t_res = scs.resample(I, int(t[-1]), t=t)
    # Estimando os parâmetros
    bounds = zip([x0[0]/10, x0[1]/10], [10*x0[0], 1000*x0[1]])
    res = differential_evolution(cost_function, list(bounds), args=(Sd_res, Id_res, y0, t_res, w))
    # Simulando os dados
    [Sa, Ia] = SIRsim(y0, t_res, res.x)
    # Save the year data
    final_data["data"]["original"].append(
        { "I": I, "B": B, "S": S, "t": t/365 + year_ref })
    final_data["data"]["resampled"].append(
        {"I": Id_res, "S": Sd_res, "t": t_res/365 + year_ref})
    final_data["data"]["simulated"].append(
        {"I": Ia, "S": Sa, "t": t_res/365 + year_ref})
    final_data["pars"]["beta"].append( res.x[0] )
    final_data["pars"]["r"].append( res.x[1] )
[22]:

r = final_data["pars"]["r"]
beta = final_data["pars"]["beta"]
years = [int(t[0]) for t in windowed_data["t"]]



p = figure(tools=TOOLS+",hover",
           y_range=(min(beta), max(beta)),
           plot_width=600, plot_height=300)


p.line(years, beta,
       legend_label="beta", line_width=4, color="#c2185b", line_cap='round', line_alpha=0.9)

p.extra_y_ranges = {"r_axis": Range1d(start=min(r), end=max(r))}
p.add_layout(LinearAxis(y_range_name="r_axis"), 'left')

p.line(years, r, y_range_name="r_axis", line_dash='dashed',
       legend_label="r", line_width=3, color="#8e44ad", line_cap='round', line_alpha=0.9)


p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.xaxis.axis_label = "Ano"
p.toolbar.autohide = True



p1 = figure(tools=TOOLS+",hover",
            x_range=p.x_range,
            plot_width=600, plot_height=300)


p1.line(data["time"], data["cases"],
        legend_label="Casos", line_width=2, color="#f4511e", line_cap='round', line_alpha=0.9)

for dataset in final_data["data"]["original"]:
    p1.line(dataset["t"], dataset["I"],
           legend_label="Casos", line_width=4, color="#f4511e", line_cap='round', line_alpha=0.9)

for dataset in final_data["data"]["simulated"]:
    p1.line(dataset["t"], dataset["I"], line_dash='dashed',
           legend_label="Estimado", line_width=3, color="#0288d1", line_cap='round', line_alpha=0.9)


p1.grid.grid_line_alpha = 0
p1.ygrid.band_fill_color = "olive"
p1.ygrid.band_fill_alpha = 0.1
p1.yaxis.axis_label = "Indivíduos"
p1.xaxis.axis_label = "Ano"
p1.toolbar.autohide = True



show(column(p,p1))

[23]:

source = ColumnDataSource(data=dict(date=data["time"].to_numpy(), value=data["cases"].to_numpy()))

p = figure(plot_height=300, plot_width=600, tools="xpan", toolbar_location=None,
           x_axis_location="above", x_range=(data["time"][5], data["time"][105]))

p.line('date', 'value', source=source,
       legend_label="Casos", line_width=2, color="#f4511e", line_cap='round', line_alpha=0.9)

for dataset in final_data["data"]["resampled"]:
    rsource = ColumnDataSource(data=dict(date=dataset["t"], value=dataset["I"]))
    p.scatter('date', 'value', source=rsource,
       legend_label="Casos reamostrados", fill_color="#8e44ad", size=4, line_alpha=0)

p.yaxis.axis_label = 'Indivíduos'
p.grid.grid_line_alpha = 0
p.xgrid.band_fill_color = "navy"
p.xgrid.band_fill_alpha = 0.1
p.toolbar.autohide = True

select = figure(title="Drag the middle and edges of the selection box to change the range above",
                plot_height=130, plot_width=600, y_range=p.y_range,
                y_axis_type=None, tools="", toolbar_location=None)

range_tool = RangeTool(x_range=p.x_range)
range_tool.overlay.fill_color = "navy"
range_tool.overlay.fill_alpha = 0.2

select.line('date', 'value', source=source,
            line_width=2, color="#f4511e", line_cap='round', line_alpha=0.9)
select.ygrid.grid_line_color = None
select.add_tools(range_tool)
select.toolbar.active_multi = range_tool
select.xaxis.axis_label = 'Ano'

show(column(p, select))

[24]:

p1 = figure(plot_height=300, plot_width=600, tools="xpan", toolbar_location=None,
            x_axis_location="above", x_range=p.x_range)

p1.line('date', 'value', source=source,
       legend_label="Casos", line_width=2, color="#f4511e", line_cap='round', line_alpha=0.9)

for dataset in final_data["data"]["simulated"]:
    psource = ColumnDataSource(data=dict(date=dataset["t"], value=dataset["I"]))
    p1.line('date', 'value', source=psource, line_dash="dashed",
           legend_label="Modelo estimado", color="#0288d1",
           line_width=4, line_cap='round', line_alpha=0.9)

    window_bound_lower = [dataset["t"][0], dataset["t"][0]]
    window_bound_upper = [dataset["t"][-1], dataset["t"][-1]]
    window_limits = [0, 10000]

    p.line(window_bound_lower, window_limits, line_dash="dashed",
           legend_label="Janela da epidemia", color="#455a64",
           line_width=2, line_cap='round', line_alpha=0.2)
    p.line(window_bound_upper, window_limits, line_dash="dashed",
           legend_label="Janela da epidemia", color="#455a64",
           line_width=2, line_cap='round', line_alpha=0.2)
    p1.line(window_bound_lower, window_limits, line_dash="dashed",
           legend_label="Janela da epidemia", color="#455a64",
           line_width=2, line_cap='round', line_alpha=0.2)
    p1.line(window_bound_upper, window_limits, line_dash="dashed",
           legend_label="Janela da epidemia", color="#455a64",
           line_width=2, line_cap='round', line_alpha=0.2)

p1.xaxis.axis_line_alpha = 0
p1.xaxis.major_label_text_color = None
p1.xaxis.major_tick_line_color = None
p1.xaxis.minor_tick_line_color = None
p1.yaxis.axis_label = 'Indivíduos'
p1.grid.grid_line_alpha = 0
p1.xgrid.band_fill_color = "navy"
p1.xgrid.band_fill_alpha = 0.1
p1.toolbar.autohide = True

show(column(p, p1, select))

Salvando as variáveis
[25]:

import pickle

save_variable = {
    "data" : data,
    "final" : final_data
}

with open('UK_estimated_data.pickle', 'wb') as handle:
    pickle.dump(save_variable, handle)

SIR 🇬🇧 Reino Unido Modelos Espaciais

Os dados

Lendo os dados de todas as cidades disponíveis:

[1]:

import os
import pandas as pd


path = "./PG_IMT/DadosEpidemia/UKCities/"
files = os.listdir(path)

dataframes = []
for file in files:
    data = pd.read_csv(path + file).to_dict()
    keys = [k for k in data.keys()]
    size = len(data[keys[0]])
    # Create a column for the city name
    city = [file[:-4] for n in range(size)]
    data["city"] = dict(zip(range(size), city))
    dataframes.append(pd.DataFrame(data))

data = pd.concat(dataframes, ignore_index=True)
data = data.dropna()
data.head()

[1]:
time cases births pop city
0 1944.016427 0 27.269231 28350.000000 Bedwellty
1 1944.054757 0 27.148291 28339.031079 Bedwellty
2 1944.093087 0 27.027352 28328.062157 Bedwellty
3 1944.131417 0 26.906413 28317.093236 Bedwellty
4 1944.169747 1 26.785473 28306.124314 Bedwellty

Determinando os parâmetros - SIR

Nesta seção, iremos utilizar a biblioteca models, que concentra todas as funções dos modelos SIR desenvolvidos nas seções anteriores.

Para uma única cidade

[2]:

from models import *

# Getting the data for a particular city
dataset = data.where(data["city"] == "London").dropna()

# Creating the data for trainning
N = int(dataset["pop"].mean())
B = dataset["births"].to_numpy()
I = dataset["cases"].to_numpy()
S = dataset["pop"].to_numpy() - I + B
t = dataset["time"].to_numpy()

# Creating the SIR model
model = ss.SIR(pop=N, focus=["I"])

# Fitting the model on data
fit_data = model.fit_multiple(S, I, B, t,
                              r_sens=[1000,10], beta_sens=[100,100])

# Plotting the model results
fig = model.result_summary(
      out_plot=True,
      plot_size=[600,400])


from bokeh.plotting import show

show(fig)

Loading BokehJS ...
Windows starting at:  [17, 46, 98, 124, 174, 222, 276, 331, 383, 429, 481]
Windows ending at:    [47, 72, 124, 153, 206, 260, 311, 361, 415, 468, 519]
Window start cases:   [41.0, 88.0, 113.0, 157.0, 308.0, 1521.0, 249.0, 183.0, 194.0, 256.0, 141.0]
New iter:::  1
         ├─ S(0) ─   2550525.415201801   I(0) ─   41.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.7605155319336778  ─  0.5771663611353828

New iter:::  2
         ├─ S(0) ─   2949891.618829986   I(0) ─   88.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.652527725167639  ─  0.5784638670213627

New iter:::  3
         ├─ S(0) ─   3302261.375966816   I(0) ─   113.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.5930110943411019  ─  0.5853040022414066

New iter:::  4
         ├─ S(0) ─   3360609.294403033   I(0) ─   157.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.5510788084784277  ─  0.5541593857111112

New iter:::  5
         ├─ S(0) ─   3379912.492617072   I(0) ─   308.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.4418921319691951  ─  0.437900484416566

New iter:::  6
         ├─ S(0) ─   3361115.5108282953   I(0) ─   1521.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.2746852290519879  ─  0.27916031863574764

New iter:::  7
         ├─ S(0) ─   3321231.492546763   I(0) ─   249.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.36599758110134534  ─  0.35877292971764685

New iter:::  8
         ├─ S(0) ─   3272868.3492476433   I(0) ─   183.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.4917754139555685  ─  0.4781307010467401

New iter:::  9
         ├─ S(0) ─   3226812.132892697   I(0) ─   194.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.5780416048451077  ─  0.5570092885781249

New iter:::  10
         ├─ S(0) ─   3199397.624806639   I(0) ─   256.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.3327289237627916  ─  0.315264257523117

New iter:::  11
         ├─ S(0) ─   3185602.2749261744   I(0) ─   141.0
         ├─ beta ─   1   r ─   0.1
         ├─ beta bound ─   0.01  ─  10
         ├─ r bound ─   0.001  ─  1.0
         └─ Defined at:  0.4244705235353749  ─  0.4056067055948321

[3]:

[3]:
3212255.4403247973

Estimando para todas cidades

[ ]:

from models import *

# Define the city names
cities = data["city"].unique()

par_data = dict()
for i, city in enumerate(cities):

    # Getting the data for a particular city
    dataset = data.where(data["city"] == city).dropna()

    # Creating the data for trainning
    N = dataset["pop"].mean()
    B = dataset["births"].to_numpy()
    I = dataset["cases"].to_numpy()
    S = dataset["pop"].to_numpy() - I + B
    t = dataset["time"].to_numpy()

    # Creating the SIR model
    model = ss.SIR(verbose=False)

    # Fitting the model on data
    par_data[city] = model.fit_multiple(S,I,B,t, out_type=1)

    # Save the summary on folder
    model.result_summary(
      save_results=True,
      plot_size=[700,500],
      folder_path="./estimation_summaries/",
      file_name=city+"_summary.png")

    print("Finished - ", city, " - ", i+1, " of ", len(cities))

Obtendo informações de Lat e Long

[ ]:

from geopy.geocoders import Nominatim

geolocator = Nominatim(user_agent="UK_EpidemicModels_App")

for city in cities:
    city_name = city.replace(".", " ")
    location = geolocator.geocode(city_name)
    par_data[city]["lat"] = location[-1][0]
    par_data[city]["lon"] = location[-1][1]
    print("- Located: ", location[-1], "  ", city_name)

Correlação dos modelos

Interpolando os modelos estimados

[ ]:

import numpy as np
import scipy.signal as scs
from scipy.interpolate import Rbf, InterpolatedUnivariateSpline

first_year = round(min([par_data[city]["time"][0] for city in cities]))
last_year = round(max([par_data[city]["time"][-1] for city in cities]))
points = int(last_year - first_year) - 1

data_struc = {
    "lat": [],
    "lon": [],
    "city": [],
    "year": [],
    "beta": [],
    "r": [],
}

points = 30
time = np.linspace(first_year, last_year, points)

for city in cities:

    # Selecionando os parâmetros estimados
    beta_res = par_data[city]["pars"]["beta"]
    r_res = par_data[city]["pars"]["r"]
    year = par_data[city]["time"]

    # Interpolando com reamostragem
    #beta_res, t_res = scs.resample(beta, points, t=year-year[0])
    #r_res, t_res = scs.resample(r, points, t=year-year[0])
    #t_res = t_res.astype(int)

    # Interpolando com Univariate Splines
    #beta_rbf = InterpolatedUnivariateSpline(year, beta_res)
    #r_rbf = InterpolatedUnivariateSpline(year, r_res)

    # Interpolando com Radial Basis Functions
    beta_rbf = Rbf(year, beta_res, function='gaussian', smooth=2)
    r_rbf = Rbf(year, r_res, function='gaussian', smooth=2)

    # Salvando do dicionário cada campo
    data_struc["lat"] += [par_data[city]["lat"] for k in range(points)]
    data_struc["lon"] += [par_data[city]["lon"] for k in range(points)]
    data_struc["city"] += [city for k in range(points)]
    data_struc["beta"] += beta_rbf(time).tolist()
    data_struc["r"]    += r_rbf(time).tolist()
    data_struc["year"] += time.tolist()

# Criando um data frame
df = pd.DataFrame(data_struc)

Visualizando a interpolação
[ ]:

# Creating the parameters plot
#

from bokeh.models   import Legend, ColumnDataSource, RangeTool, LinearAxis, Range1d
from bokeh.palettes import brewer, Inferno256
from bokeh.plotting import figure, show
from bokeh.layouts  import column
from bokeh.io       import output_notebook


p_beta = figure(
    title="Beta Parameter",
    y_axis_type="log",
    plot_width=700,
    plot_height=500
)

p_r = figure(
    title="R Parameter",
    y_axis_type="log",
    plot_width=700,
    plot_height=500
)

legend_it, legend_it_r = [], []
for i, city in enumerate(cities):

    color = Inferno256[int((i/len(cities))*256)]

    df_filt = df.where(df["city"]==city).dropna()

    c = p_beta.line(
        df_filt["year"].to_list(),
        abs(df_filt["beta"]),
        line_width=4,
        line_cap="round",
        color=color
    )

    cr = p_r.line(
        df_filt["year"],
        abs(df_filt["r"]),
        line_width=4,
        line_cap="round",
        color=color
    )

    legend_it.append((city, [c]))
    legend_it_r.append((city, [cr]))

p_beta.grid.grid_line_alpha = 0
p_beta.ygrid.band_fill_color = "olive"
p_beta.ygrid.band_fill_alpha = 0.1
p_beta.xaxis.axis_label = "Ano"

legend = Legend(items=legend_it, location=(0, -10))
legend.click_policy="mute"
p_beta.add_layout(legend, "right")
p_beta.legend.click_policy="hide"

p_r.grid.grid_line_alpha = 0
p_r.ygrid.band_fill_color = "olive"
p_r.ygrid.band_fill_alpha = 0.1
p_r.xaxis.axis_label = "Ano"

legend = Legend(items=legend_it_r, location=(0, -10))
legend.click_policy="mute"
p_r.add_layout(legend, "right")
p_r.legend.click_policy="hide"

show(column(p_beta, p_r))

Criando os data frames para matrix de correlação

[ ]:

r_data = {}
beta_data = {}

for city in cities:

    df_filt = df.where(df["city"]==city).dropna()

    data_dict = dict()
    zipped_data = zip(
        df_filt["year"].to_list(),
        df_filt["beta"].to_list())
    for y, b in zipped_data:
        data_dict[y] = b
    beta_data[city] = data_dict

    data_dict = dict()
    zipped_data = zip(
        df_filt["year"].to_list(),
        df_filt["r"].to_list())
    for y, r in zipped_data:
        data_dict[y] = r
    r_data[city] = data_dict

r_df = pd.DataFrame(r_data)
r_df = r_df.sort_index()
r_df.index = (1000 * r_df.index.to_numpy()).astype(int)

beta_df = pd.DataFrame(beta_data)
beta_df = beta_df.sort_index()
beta_df.index = (1000 * beta_df.index.to_numpy()).astype(int)
beta_df.head()

Correlação do parâmetro \(\beta\)

[ ]:

import seaborn as sns
import matplotlib.pyplot as plt

# Correlation Matrix Heatmap
f, ax = plt.subplots(figsize=(18, 15))

beta_corr = beta_df.corr()
hm = sns.heatmap(round(beta_corr,2),
                 annot=True,
                 ax=ax,
                 cmap="coolwarm",
                 fmt='.2f',
                 linewidths=.05)

f.subplots_adjust(top=0.93)

t = f.suptitle('Beta Parameter - Features Correlation Heatmap', fontsize=14)

Correlação do parâmetro \(r\)

[ ]:

# Correlation Matrix Heatmap
f, ax = plt.subplots(figsize=(18, 15))

r_corr = r_df.corr()
hm = sns.heatmap(round(r_corr,2),
                 annot=True,
                 ax=ax,
                 cmap="coolwarm",
                 fmt='.2f',
                 linewidths=.05)

f.subplots_adjust(top=0.93)

t = f.suptitle('r Parameter - Features Correlation Heatmap', fontsize=14)

Correlação do modelo

[ ]:

f, ax = plt.subplots(figsize=(18, 15))

model_corr = 0.5 * ( r_corr +  beta_corr )
hm = sns.heatmap(round(model_corr,2),
                 annot=True,
                 ax=ax,
                 cmap="coolwarm",
                 fmt='.2f',
                 linewidths=.05)

f.subplots_adjust(top=0.93)

t = f.suptitle('Model - Features Correlation Heatmap', fontsize=14)

SIR 🇧🇷 Parâmetros variantes

Os Dados

[1]:

import pandas as pd

# reading csv file
BRdata = pd.read_csv("./PG_IMT/DadosEpidemia/CoVid19.csv")
BRdata.head()

[1]:
date country state city newDeaths deaths newCases totalCases deathsMS totalCasesMS deaths_per_100k_inhabitants totalCases_per_100k_inhabitants deaths_by_totalCases recovered suspects tests tests_per_100k_inhabitants
0 2020-02-25 Brazil SP TOTAL 0 0 1 1 0 0 0.0 0.00218 0.0 NaN NaN NaN NaN
1 2020-02-25 Brazil TOTAL TOTAL 0 0 1 1 0 0 0.0 0.00048 0.0 NaN NaN NaN NaN
2 2020-02-26 Brazil SP TOTAL 0 0 0 1 0 1 0.0 0.00218 0.0 NaN NaN NaN NaN
3 2020-02-26 Brazil TOTAL TOTAL 0 0 0 1 0 1 0.0 0.00048 0.0 NaN NaN NaN NaN
4 2020-02-27 Brazil SP TOTAL 0 0 0 1 0 1 0.0 0.00218 0.0 NaN NaN NaN NaN

Filtrando e limpando os dados

[2]:

is_SP =  BRdata['state']=='TOTAL'
SPdata = BRdata[is_SP]

# SPdata = SPdata.fillna(0)

SPLimpo = SPdata[SPdata.recovered.notnull()]
SPLimpo.head()

[2]:
date country state city newDeaths deaths newCases totalCases deathsMS totalCasesMS deaths_per_100k_inhabitants totalCases_per_100k_inhabitants deaths_by_totalCases recovered suspects tests tests_per_100k_inhabitants
321 2020-03-23 Brazil TOTAL TOTAL 9 34 358 1952 34 1891 0.01618 0.92887 0.01742 8.0 15867.0 NaN NaN
349 2020-03-24 Brazil TOTAL TOTAL 13 47 303 2255 46 2201 0.02237 1.07306 0.02084 20.0 17700.0 NaN NaN
377 2020-03-25 Brazil TOTAL TOTAL 12 59 311 2566 57 2433 0.02808 1.22105 0.02299 27.0 27227.0 NaN NaN
405 2020-03-26 Brazil TOTAL TOTAL 18 77 424 2990 77 2915 0.03664 1.42281 0.02575 42.0 48793.0 NaN NaN
433 2020-03-27 Brazil TOTAL TOTAL 16 93 486 3476 92 3417 0.04425 1.65408 0.02675 42.0 50684.0 NaN NaN
[3]:

from datetime import datetime

first_date = SPLimpo["date"].iloc[0]
first_date = datetime.fromisoformat(first_date)

SPsir = SPLimpo[["totalCases", "deaths", "recovered"]].to_numpy()

SPsir[:,0:]
SPI = SPsir[:,0]
auxM = SPsir[:,1]
SPM = auxM
auxR = SPsir[:,2]
SPR = auxR

Visualizando a evolução

[4]:

from bokeh.models   import Legend, ColumnDataSource, RangeTool, LinearAxis, Range1d
from bokeh.palettes import brewer, Inferno256
from bokeh.plotting import figure, show
from bokeh.layouts  import column
from bokeh.io       import output_notebook

output_notebook()

from datetime import timedelta

# Criando o vetor de tempo
t_num = range(len(SPM))
date_vec = [ first_date + timedelta(days=k) for k in t_num]


# Criando a figura
p = figure(
    tools="hover",
    title="Evolução do COVID",
    x_axis_type='datetime',
    y_axis_type="log",
    plot_width=650,
    plot_height=500
)

# Incluindo os dados de mortes
p.line(
    date_vec, SPM,
    legend_label="Mortes",
    line_width=4,
    line_cap="round",
    color="#de425b"
)

# Incluindo os dados de infectados
p.line(
    date_vec, SPI,
    legend_label="Infectados",
    line_width=4,
    line_cap="round",
    color="#ffd885"
)

# Incluindo os dados de recuperados
p.line(
    date_vec, SPR,
    legend_label="Removidos",
    line_width=4,
    line_cap="round",
    color="#99d594"
)


p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"
p.legend.click_policy="hide"
p.legend.location = "bottom_right"

show(p)

Loading BokehJS ...

Análise de um modelo com \(\beta\) variável no tempo

Uma outra forma de abordar a modelagme de dados é assumindo que a taxa de transmissão \(\beta\) é variante no tempo. Vamos, neste sentido, utilizar a abordagem elaborada por Mark Pollicott, Hao Wang & Howard (Howie) Weiss em 2012 - Extracting the time-dependent transmission rate from infection data via solution of an inverse ODE problem

Reescrevendo o conjunto de equações diferenciais que caracteriza o modelo é descrito abaixo. No modelo $:nbsphinx-math:beta`(t) - :nbsphinx-math:text{representa a taxa de transmissão ou taxa efetiva de contato}` $ e \(r - \text{a taxa de remoção ou recuperação.}\)

\[\begin{split}\begin{split} \frac{dS(t)}{dt} & = -\mathbf{\beta(t)} S(t) I(t) \\ \frac{dI(t)}{dt} & = \mathbf{\beta(t)} S(t) I(t) - rI(t) \\ \frac{dR(t)}{dt} & = r I(t) \end{split}\end{split}\]

Gostaríamos de identificar quais parâmetros \(\beta(t)\) e \(r\) resultam num melhor ajuste do modelo para os dados de S,I e R

Proposta de algoritmo para extração do \(\beta(t)\) a partir dos dados

  1. Interpolar os dados para suvizar eventuais variações gerando uma função \(f(t)\)
  2. Verificar a seguinte condição: \(\frac{f'(t)}{f(t)}>-r\). Neste algoritmo a taza de romação \(r\) é considerada constante.
  3. Determinando a função \(p(t) = \frac{f''(t)f(t)-f'^2(t)}{f(t)(f'(t)+rf(t))}\).
  4. Escolhendo um valor inicial para \(\beta(0)\) e determinando a integral \(P(t) = \int_0^t p(\tau) d\tau\).
  5. Verificar a seguinte condição: \(\beta(0) <\frac{1}{\int_0^t e^{P(s)}f(s) ds}\).
  6. O parâmetro \(\mathbf{\beta(t)}\) pode ser calculado por: \(\beta(t) = \frac{1}{\frac{e^{-P(t)}}{\beta(0)}-e^{-P(t)}\int_0^t e^{P(s)}f(s) ds}\).
1) Interpolar os dados para suvizar eventuais variações gerando uma função \(f(t)\)
[5]:

import numpy as np
from scipy.interpolate import Rbf, UnivariateSpline

# Generating weights for polynomial function with degree =2
weights = np.polyfit(t_num, SPI, 2)
rbf = Rbf(t_num, SPI, function='gaussian')
spline = UnivariateSpline(t_num, SPI)
print(weights)


# Generating model with the given weights
model = np.poly1d(weights)

t_sim = np.linspace(t_num[0], t_num[-1], 100)

# Prediction on validation set
pred = model(t_sim)
pred_rbf = rbf(t_sim)
pred_spline = spline(t_sim)

# We will plot the graph for 70 observations only

# Criando a figura
p = figure(
    tools="hover",
    title="Dados interpolados",
    x_axis_type='datetime',
#     y_axis_type="log",
    plot_width=650,
    plot_height=500
)

# Incluindo os dados de mortes
p.line(
    t_num, SPI,
    legend_label="Medidas",
    line_width=4,
    line_cap="round",
    color="#de425b"
)

# Incluindo os dados de infectados
p.line(
    t_sim, pred,
    legend_label="Interpolação",
    line_width=4,
    line_cap="round",
    color="#1e88e5"
)

p.line(
    t_sim, pred_rbf,
    legend_label="Interpolação - RBF",
    line_width=3,
    line_dash="dashed",
    line_cap="round",
    color="#512da8"
)

p.line(
    t_sim, pred_spline,
    legend_label="Interpolação - Spline",
    line_width=3,
    line_dash="dashed",
    line_cap="round",
    color="#00796b"
)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"
p.legend.click_policy="hide"
p.legend.location = "bottom_right"

show(p)


[  67.47775863 -501.31086778 4785.98379386]

SIR 🇧🇷 Brasil - Dados Reais

Os dados

Lendo e visualizandos os dados do Brasil do COVID fornecidos por

Os dados são atualizados diariamente… Então normalmente a cada dia esse pipeline pode mudar seus resultados.

[1]:

import pandas as pd

# Lendo os dados online - wcota
data_path = 'https://raw.githubusercontent.com/wcota/covid19br/master/cases-brazil-states.csv'

online_data = pd.read_csv(data_path, delimiter=",")
online_data.head()

[1]:
date country state city newDeaths deaths newCases totalCases deathsMS totalCasesMS deaths_per_100k_inhabitants totalCases_per_100k_inhabitants deaths_by_totalCases recovered suspects tests tests_per_100k_inhabitants
0 2020-02-25 Brazil SP TOTAL 0 0 1 1 0 0 0.0 0.00218 0.0 NaN NaN NaN NaN
1 2020-02-25 Brazil TOTAL TOTAL 0 0 1 1 0 0 0.0 0.00048 0.0 NaN NaN NaN NaN
2 2020-02-26 Brazil SP TOTAL 0 0 0 1 0 1 0.0 0.00218 0.0 NaN NaN NaN NaN
3 2020-02-26 Brazil TOTAL TOTAL 0 0 0 1 0 1 0.0 0.00048 0.0 NaN NaN NaN NaN
4 2020-02-27 Brazil SP TOTAL 0 0 0 1 0 1 0.0 0.00218 0.0 NaN NaN NaN NaN

Filtrando e limpando os dados

[2]:

selected_state = "SP"

at_state = online_data['state']==selected_state
local_data = online_data[at_state]
local_data = local_data[local_data.recovered.notnull()]
#local_data = local_data.fillna(method="backfill")

local_data.head()

[2]:
date country state city newDeaths deaths newCases totalCases deathsMS totalCasesMS deaths_per_100k_inhabitants totalCases_per_100k_inhabitants deaths_by_totalCases recovered suspects tests tests_per_100k_inhabitants
347 2020-03-24 Brazil SP TOTAL 10 40 65 810 40 810 0.08711 1.76397 0.04938 1.0 4572.0 NaN NaN
375 2020-03-25 Brazil SP TOTAL 8 48 52 862 48 862 0.10453 1.87722 0.05568 1.0 4300.0 NaN NaN
403 2020-03-26 Brazil SP TOTAL 10 58 191 1053 58 1052 0.12631 2.29317 0.05508 1.0 14312.0 NaN NaN
431 2020-03-27 Brazil SP TOTAL 10 68 170 1223 68 1223 0.14809 2.66338 0.05560 1.0 14312.0 NaN NaN
459 2020-03-28 Brazil SP TOTAL 16 84 183 1406 84 1406 0.18293 3.06191 0.05974 1.0 14312.0 NaN NaN
[3]:

import numpy as np
from datetime import datetime

first_date = local_data["date"].iloc[0]
first_date = datetime.fromisoformat(first_date)

if selected_state == "SP":
#     N = 11869660
    N = 44.01e6
elif selected_state == "TOTAL":
    N = 220e6

I = list()                                       # <- I(t)
R = local_data["recovered"].iloc[1:].to_numpy()  # <- R(t)
M = local_data["newDeaths"].iloc[1:].to_numpy()  # <- M(t)
nR = np.diff(local_data["recovered"].to_numpy()) # <- dR(t)/dt
nC = local_data["newCases"].iloc[1:].to_numpy()  # <- nC(t)/dt

I = [ local_data["totalCases"].iloc[1] ]         # I(0)

# I(t) <- I(t-1) + newCases(t) - newMortes(t) - newRecovered(t)
for t in range(len(M)-1):
    I.append(I[-1] + nC[t] - M[t] - nR[t])
I = np.array(I)

Visualizando a evolução

[18]:

from bokeh.models   import Legend, ColumnDataSource, RangeTool, LinearAxis, Range1d, HoverTool
from bokeh.palettes import brewer, Inferno256
from bokeh.plotting import figure, show
from bokeh.layouts  import column
from bokeh.io       import output_notebook

output_notebook()

from datetime import timedelta

# Criando o vetor de tempo
date_vec = [ first_date + timedelta(days=k) for k in range(len(M))]

# Criando os valores para legenda no plot
year =  [str(int(d.year)) for d in date_vec ]
month = [("0"+str(int(d.month)))[-2:] for d in date_vec ]
day =   [("0"+str(int(d.day)))[-2:] for d in date_vec ]



# Criando a fonte de dados
source = ColumnDataSource(data={
    'Data'       : date_vec,
    'd': day, 'm': month, 'y': year,
    'Infectados' : I,
    'Removidos'  : R,
    'Mortes'     : M,
})


# Criando a figura
p = figure(plot_height=500,
           plot_width=600,
           x_axis_type="datetime",
           tools="",
           toolbar_location=None,
           y_axis_type="log",
           title="Evolução do COVID - São Paulo")

# Preparando o estilo
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

# Incluindo as curvas
i_p = p.line(x='Data', y='Infectados', legend_label="Infectados", line_cap="round", line_width=3, color="#ffd885", source=source)
m_p = p.line(x='Data', y='Mortes', legend_label="Mortes", line_cap="round", line_width=3, color="#de425b", source=source)
r_p = p.line(x='Data', y='Removidos', legend_label="Removidos", line_cap="round", line_width=3, color="#99d594", source=source)

# Colocando as legendas
p.legend.click_policy="hide"
p.legend.location = "top_left"

# Incluindo a ferramenta de hover
p.add_tools(HoverTool(
    tooltips=[
        ( 'Indivíduos', '$y{i}'),
        ( 'Data',       '@d/@m/@y' ),
    ],
    renderers=[
        r_p, i_p, m_p
    ]
))

show(p)

Loading BokehJS ...

O problema

O conjunto de equações diferenciais que caracteriza o modelo é descrito abaixo. No modelo $:nbsphinx-math:beta `- :nbsphinx-math:text{representa a taxa de transmissão ou taxa efetiva de contato}` $ e \(r - \text{a taxa de remoção ou recuperação.}\)

\[\begin{split}\begin{split} \frac{dS(t)}{dt} & = -\beta S(t) I(t) \\ \frac{dI(t)}{dt} & = \beta S(t) I(t) - rI(t) \\ \frac{dR(t)}{dt} & = r I(t) \end{split}\end{split}\]

Gostaríamos de identificar quais parâmetros \(\beta\) e \(r\) resultam num melhor ajuste do modelo para os dados de S,I e R

[5]:

# Importando o modelo SIR
from models import *

sir_model = ss.SIR(pop=N, focus=["I", "R"])

Loading BokehJS ...

Estimando os parâmetros

Para estimarmos os parâmetros do modelo \(\mathbf{\beta}\) e \(\mathbf{r}\), vamos utilizar inicialmente o método de mínimos quadrados. Podemos então formular o problema a partir da Equação abaixo. Na Equação \(y_m(k)\) representa o dado real em cada amostra \(k\); \(y_s(\theta,k)\) representa o valor estimado a partir da simulação do modelo para uma determinada amostra \(k\) e \(\theta\) representa o vetor ed parâmetros \(\theta = [ \beta \; \; r]^T\).

\[min_{\theta}= \sum_{k=1}^{K}(y_m(k) - y_s(\theta,k))^2\]

A equação formula a pergunta: quais os valores de \(beta\) e \(r\) que minizam o erro quadrático quando comparados com os dados reais.

[12]:

import numpy as np

S = N - I - R

time = np.linspace(0, len(I), len(I))

# Estimando os parâmetros
sir_model.fit(S, I, R, time, sample_ponder=True, resample=True, beta_sens=[1000,10], r_sens=[1000,10])

r_included = True

         ├─ Resample from sizes ─  65 65 65 65
         └─ Resample to sizes ─    1540 1540 1540 1540
         ├─ S(0) ─ I(0) ─ R(0) ─  [44009136.41827327, 862.5810526571099, 1.0006740750082481]
         ├─ beta ─   1   r ─   0.14285714285714285
         ├─ beta bound ─   0.001  ─  10
         ├─ r bound ─   0.00014285714285714284  ─  1.4285714285714284
         ├─ equation weights ─   [0.0016004861866103122, 1, 1]
         └─ Defined at:  0.08388235287485474  ─  0.010752489662406836

[13]:
sir_model.parameters[0]/sir_model.parameters[1] # Ro <- \beta * (1 / \r)
[13]:
7.801202838457649
[19]:

# Criando a figura
p1 = figure(plot_height=500,
           plot_width=600,
           x_axis_type="datetime",
           tools="",
           toolbar_location=None,
#            y_axis_type="log",
           title="Evolução do COVID - São Paulo")

# Preparando o estilo
p1.grid.grid_line_alpha = 0
p1.ygrid.band_fill_color = "olive"
p1.ygrid.band_fill_alpha = 0.1
p1.yaxis.axis_label = "Indivíduos"
p1.xaxis.axis_label = "Dias"

# Incluindo as curvas
p1.line(time, I, legend_label="Infectados", line_cap="round", line_width=3, color="#ffd885")
#p1.line(time, , legend_label="Mortes", line_cap="round", line_width=3, color="#de425b")
p1.line(time, R, legend_label="Removidos", line_cap="round", line_width=3, color="#99d594")

p1.scatter(sir_model.pipeline["resample"]["after"]["t"],
           sir_model.pipeline["resample"]["after"]["I"],
           marker="circle",line_color="#6666ee", fill_color="#ee6666", fill_alpha=0.5, size=3)

p1.scatter(sir_model.pipeline["resample"]["after"]["t"],
           sir_model.pipeline["resample"]["after"]["R"],
           marker="circle",line_color="#6666ee", fill_color="#ee6666", fill_alpha=0.5, size=3)

# Colocando as legendas
p1.legend.click_policy="hide"
p1.legend.location = "top_left"


show(p1)


[20]:

if r_included:
    initial = [S[0], I[0], R[0]]
else:
    initial = [S[0], I[0]]

results = sir_model.predict(initial, time)

[21]:

# Incluindo os dados de infectados
im_p = p.line(
    date_vec, results[1],
    legend_label="Infectados - Modelo",
    line_width=4,
    line_dash="dashed",
    line_cap="round",
    color="#f57f17"
)

# Incluindo os dados de recuperados
if r_included:
    rm_p = p.line(
        date_vec, results[2],
        legend_label="Removidos - Modelo",
        line_dash="dashed",
        line_width=4,
        line_cap="round",
        color="#1b5e20"
    )

show(p)

Predições utilizando o modelo

[17]:

# Criando os valores de tempo para previsão - 120 dias
t_sim = np.linspace(0, len(I) + 120, len(I) + 120)
date_vec_sim = [first_date + timedelta(days=k) for k in t_sim]

# Prevendo para os valores selecionados
prediction = sir_model.predict(initial, t_sim)



# Criando o gráfico com as predições

# Criando os valores para legenda no plot
year_sim =  [str(int(d.year)) for d in date_vec_sim ]
month_sim = [("0"+str(int(d.month)))[-2:] for d in date_vec_sim ]
day_sim =   [("0"+str(int(d.day)))[-2:] for d in date_vec_sim ]

if r_included:
    accum_Infect = [0]
    for i in N - prediction[1] - prediction[2]:
        accum_Infect.append(accum_Infect[-1]+i)

    accum_Infect = sir_model.parameters[1] * np.array(accum_Infect)

# Criando a fonte de dados
if r_included:
    source = ColumnDataSource(data={
        'Data'       : date_vec,
        'd': day, 'm': month, 'y': year,
        'Infectados' : I,
        'Removidos'  : R,
        'Mortes'     : M,
        'InfecModelo' : prediction[1],
        'RemovModelo' : prediction[2],
        'AccumInfect' : accum_Infect,
        'SucetModelo' : N - prediction[1] - prediction[2],
        'DataModelo'  : date_vec_sim,
        'ds': day_sim, 'ms': month_sim, 'ys': year_sim
    })
else:
    source = ColumnDataSource(data={
        'Data'       : date_vec,
        'd': day, 'm': month, 'y': year,
        'Infectados' : I,
        'Removidos'  : R,
        'Mortes'     : M,
        'InfecModelo' : prediction[1],
        'DataModelo'  : date_vec_sim,
        'ds': day_sim, 'ms': month_sim, 'ys': year_sim
    })


# Criando a figura
p = figure(plot_height=700,
           plot_width=800,
           x_axis_type="datetime",
           tools="",
           toolbar_location=None,
           y_axis_type="log",
           title="Previsão do COVID - Brasil")

# Preparando o estilo
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

# Incluindo as curvas
i_p = p.line(x='Data', y='Infectados', legend_label="Infectados", line_cap="round", line_width=3, color="#ffd885", source=source)
m_p = p.line(x='Data', y='Mortes', legend_label="Mortes", line_cap="round", line_width=3, color="#de425b", source=source)
r_p = p.line(x='Data', y='Removidos', legend_label="Removidos", line_cap="round", line_width=3, color="#99d594", source=source)

mp_p = p.line(x='DataModelo', y='InfecModelo', legend_label="Infectados - Modelo", line_dash="dashed", line_cap="round", line_width=4, color="#f57f17", source=source)

renders = [i_p, m_p, r_p, mp_p]

if r_included:
    rp_p = p.line(x='DataModelo', y='RemovModelo', legend_label="Removidos - Modelo", line_dash="dashed", line_cap="round", line_width=4, color="#1b5e20", source=source)
    renders.append(rp_p)

# Colocando as legendas
p.legend.click_policy="hide"
p.legend.location = "top_left"

# Incluindo a ferramenta de hover
p.add_tools(HoverTool(
    tooltips=[
        ( 'Indivíduos', '$y{0.00 a}' ),
        ( 'Data',       '@ds/@ms/@ys'),
    ],
    renderers=renders
))

show(p)

BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('AccumInfect', 186), ('Data', 65), ('DataModelo', 185), ('InfecModelo', 185), ('Infectados', 65), ('Mortes', 65), ('RemovModelo', 185), ('Removidos', 65), ('SucetModelo', 185), ('d', 65), ('ds', 185), ('m', 65), ('ms', 185), ('y', 65), ('ys', 185)

SIR 🇨🇳 China - Dados Reais

Os dados

Para coletar os dados de outros países, foram utilizadas as APIs dos sites:

[1]:

import requests
import pandas as pd

covid_api = 'https://corona-api.com/countries/'
rest_countries = 'https://restcountries.eu/rest/v2/alpha/'
country = 'CN' # Alpha-2 ISO3166

data_json =  requests.get(covid_api + country).json()
country = requests.get(covid_api + country).json()

N = country['data']['population']

print(country['data']['name'])

China

Organizando os dados

[2]:

from datetime import datetime

df = pd.DataFrame(data_json['data']['timeline'])
df = df.sort_values('date').reset_index()

from datetime import datetime, timedelta
df['date'] = [datetime.fromisoformat(f) for f in df['date']]
df = df.drop_duplicates(subset='date', keep = 'last')

# Criando o vetor de tempo
first_date = df['date'].iloc[0]
size_days = (df['date'].iloc[-1] - df['date'].iloc[0]).days
date_vec = [first_date + timedelta(days=k) for k in range(size_days)]

new_df = pd.DataFrame(date_vec, columns=['date'])
new_df = pd.merge(new_df, df, how='left', on= 'date')
new_df = new_df.drop(columns= ['index',  'updated_at', 'is_in_progress'])

for col in new_df.columns[1:]:
    new_df[col] = new_df[col].interpolate(method='polynomial', order=1)
df = new_df.dropna()

df.head()

[2]:
date deaths confirmed active recovered new_confirmed new_recovered new_deaths
0 2020-01-21 17 548 503 28 548 28 17
1 2020-01-22 18 643 595 30 95 2 1
2 2020-01-23 26 920 858 36 277 6 8
3 2020-01-24 42 1406 1325 39 486 3 16
4 2020-01-25 56 2075 1970 49 669 10 14

Visualizando os dados

[3]:

from bokeh.models   import Legend, ColumnDataSource, RangeTool, LinearAxis, Range1d, HoverTool
from bokeh.palettes import brewer, Inferno256
from bokeh.plotting import figure, show
from bokeh.layouts  import column
from bokeh.io       import output_notebook

output_notebook()

import numpy as np

# Criando os valores para legenda no plot
year =  [str(int(d.year)) for d in df['date'] ]
month = [("0"+str(int(d.month)))[-2:] for d in df['date'] ]
day =   [("0"+str(int(d.day)))[-2:] for d in df['date'] ]

# Criando a fonte de dados
source = ColumnDataSource(data={
    'Data'       : df['date'].values,
    'd': day, 'm': month, 'y': year,
    'Infectados Acc' : df['confirmed'].values,
    'Mortes'     : df['deaths'].values,
    'Ativo'     : df['active'].values,
    'Recuperados': df['recovered'].values
})


# Criando a figura
p = figure(plot_height=500,
           plot_width=600,
           x_axis_type="datetime",
           tools="",
           #y_axis_type="log",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

# Incluindo as curvas
i_p = p.line(x='Data', y='Ativo',
             legend_label="Infectados",
             line_cap="round", line_width=5, color="#c62828", source=source)
m_p = p.line(x='Data', y='Mortes',
             legend_label="Mortes",
             line_cap="round", line_width=5, color="#512da8", source=source)
c_p = p.line(x='Data', y='Infectados Acc',
             legend_label="Infectados Acc",
             line_cap="round", line_width=5, color="#0288d1", source=source)
r_p = p.line(x='Data', y='Recuperados',
             legend_label="Recuperados",
             line_cap="round", line_width=5, color="#388e3c", source=source)

# Colocando as legendas
p.legend.click_policy="hide"
# p.legend.location = "top_left"
p.legend.location = "top_left"

# Incluindo a ferramenta de hover
p.add_tools(HoverTool(
    tooltips=[
        ( 'Indivíduos', '$y{i}'),
        ( 'Data',       '@d/@m/@y' ),
    ],
    renderers=[
        m_p, i_p, c_p, r_p
    ]
))

show(p)

Loading BokehJS ...
Verificando os dados
[4]:

dif_I = np.diff(df['active'])
cum = []
cum.append(dif_I[0])

for k, i in enumerate(dif_I):
    cum.append(cum[-1] + i)

cum = np.array(cum)
cum += df['deaths'].to_numpy() + df['recovered'].to_numpy()

print("Erro entre casos acumulados e valores de confirmados: {}".format(
        round(sum((cum - df['confirmed'].values)**2 / len(cum)),2) ) )

Erro entre casos acumulados e valores de confirmados: 168921.0
Criando os dados SIR
[5]:

I = df['active'].to_numpy()
R = df['recovered'].to_numpy()
M = df['deaths'].to_numpy()
S = N - R - I

# Creating the time vector
t = np.linspace(0, len(I), len(I))


Sd, Id, Md, Rd, td = S, I, M, R, t

[7]:
Ro = (np.log(S[0]/N) - np.log(S[-1]/N))/(S[0]/N - S[-1]/N)
print("Calculo do Ro baseado nos dados:", Ro)
Calculo do Ro baseado nos dados: 1.0000300994523579

Estimando utilizando todos os dados

[8]:

from models import *

dataset = dict(S=Sd, I=Id, R=Rd)

# Create the model
sir_model = ss.SIR(pop=N, focus=["S", "I", "R"])

# Adjust the parameters
sir_model.fit(dataset, td,
              search_pop=True,
              pop_sens=[0.00001,0.05],
              beta_sens=[100,10000],
              r_sens=[100,100])


# Predict the model
sim_res = sir_model.predict((Sd[0],Id[0],Rd[0]), td)

Loading BokehJS ...
         ├─ S(0) ─ I(0) ─ R(0) ─  [1330043469, 503, 28]
         ├─ beta ─   1   r ─   0.14285714285714285
         ├─ beta bound ─   0.01  ─  10000
         ├─ r bound ─   0.0014285714285714286  ─  14.285714285714285
         ├─ equation weights ─   [7.51893117169047e-10, 7.90340642628301e-05, 1.8154828650455473e-05]
         ├─ Running on ─  differential_evolution SciPy Search Algorithm
         └─ Defined at:  5100.406034312756  ─  0.04454555270373974

[9]:

print("Parâmetros estimados: ", sir_model.parameters)
print("Suposto Ro: ", sir_model.parameters[0] * sir_model.parameters[-1] / sir_model.parameters[1])
print("Dias contaminados: ", 1 / sir_model.parameters[1])

Parâmetros estimados:  [5.10040603e+03 4.45455527e-02 6.14351607e-05]
Suposto Ro:  7.034243490574524
Dias contaminados:  22.448930124421754
[14]:

p = figure(plot_height=500,
           plot_width=800,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

p.line(t, I,
       legend_label="Infectados", color="#ff6659", line_width=4)
p.line(t, R,
       legend_label="Removidos", color="#76d275", line_width=4)

# Show the results
p.line(td, sim_res[1],
       legend_label="Infectados - Modelo", line_dash="dashed", color="#d32f2f", line_width=3)
p.line(td, sim_res[2],
    legend_label="Removidos - Modelo", line_dash="dashed", color="#43a047", line_width=3)
p.line(td, sim_res[0],
      legend_label="Suscetíveis Ponderados - Modelo", line_dash="dashed", color="#1e88e5", line_width=3)
p.add_layout(p.legend[0], 'right')

show(p)

Monte Carlo

Nesta parte, faremos um teste aumentando a quantidade de amostras de treinamento e prevendo o momento do pico da epidemia a medida que mais dias são utilizados para treinamento. Esse estudo vai possibilitar a análise da certeza da previsão do pico da epidemia antes desse acontecer.

[15]:

saved_param = {'r':[], 'beta':[], 'pop':[]}
saved_prediction = []

start_day = 10
pred_t = np.array(range(int(td[-1])))

for i in range(start_day, len(I)):

    dataset = dict(S=Sd[:i], I=Id[:i], R=Rd[:i])
    td_ = td[:i]

    # Create the model
    sir_model = ss.SIR(pop=N, focus=["S", "I", "R"], verbose = False)

    # Adjust the parameters
    sir_model.fit(dataset, td_,
                  search_pop=True,
                  pop_sens=[0.00001,0.05],
                  beta_sens=[100,10000],
                  r_sens=[100,100])

    saved_param['beta'].append(sir_model.parameters[0])
    saved_param['r'].append(sir_model.parameters[1])
    saved_param['pop'].append(sir_model.parameters[2])

    saved_prediction.append(sir_model.predict((Sd[0],Id[0], Rd[0]), pred_t))
/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
/opt/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py:248: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
[17]:
import pickle

with open("./China_II_mc_runs.pickle", "wb") as handle:
    pickle.dump({"pars":saved_param, "pred":saved_prediction}, handle)

Análise do uso do sistema de saúde

Nesta análise, mostramos o erro percentual do quanto antes do pico, conseguimos prever a quantidade de pessoas que realmente serão identificadas como infectadas, uma vez que \(R(\infty)\) é a quantidade de pessoas recuperadas totais, daquelas que foram noficadas como infectadas no sistema de saúde. Desta forma segue o erro proporcional do erro a medida em que novos dados diários foram incluidos no modelo:

[34]:

x = range(start_day, len(I))
usage_error = [ 100 * abs(p*N - Rd[-1]) / Rd[-1] for p in saved_param['pop']]

fig9 = go.Figure()

fig9.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=usage_error[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="ε(%) = %{y:.0f}, <br> com %{x:.0f} dias de dados."))
fig9.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=usage_error[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="ε(%) = %{y:.0f}, <br> com %{x:.0f} dias de dados."))
fig9.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[usage_error[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um ε(%) = %{y:.0f}."))
fig9.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[usage_error[-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para ε(%) = %{y:.4f}."))
fig9.update_layout(template='xgridoff', yaxis_range=[-1,100],
                  legend_orientation="h", legend=dict(x=0.10, y=1.05),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Erro (%) do R(∞) - notificados',
                  title_text="Erro da previsão do uso do sistema de saúde - 🇨🇳 " + country['data']['name'])

fig9.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_China_Model_Analyser_19_0.png

Visualizando as previsões de \(I(t)\)

Vamos analisar as previsões quando somente os dados antes do pico são fornecidos ao modelo, e as previões utilizando os dados após o pico:

[19]:

visual_peak = 30

p2 = figure(plot_height=500,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p2.grid.grid_line_alpha = 0
p2.ygrid.band_fill_color = "olive"
p2.ygrid.band_fill_alpha = 0.1
p2.yaxis.axis_label = "Indivíduos"
p2.xaxis.axis_label = "Dias"

# Incluindo as curvas
p2.line(td, Id,
       legend_label="Infectados",
       line_cap="round", line_width=3, color="#c62828")

for data in saved_prediction[visual_peak-start_day:]:
    p2.line(pred_t[:len(td)], -data[1][:len(td)],
           legend_label="Previsão Infectados - Depois do pico",
           line_cap="round", line_dash="dashed", line_width=4, color="#ffa000", line_alpha = 0.1)

for data in saved_prediction[:visual_peak-start_day]:
    p2.line(pred_t[:len(td)], -data[1][:len(td)],
           legend_label="Previsão Infectados - Antes do pico",
           line_cap="round", line_dash="dashed", line_width=4, color="#42a5f5", line_alpha = 0.3)

# Colocando as legendas
p2.legend.click_policy="hide"
p2.legend.location = "top_right"

show(p2)

Visualizando os ajustes dos grupos

Aqui vamos analisar as previsões obtidas para cada grupo \(S(t)\), \(I(t)\) e \(R(t)\), a medida que mais dias foram informados ao modelo.

[26]:

p3 = figure(plot_height=350,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - I(t) e R(t) - " + country['data']['name'])

p4 = figure(plot_height=350,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - S(t) - " + country['data']['name'])

plot_all = True

# Preparando o estilo
p3.grid.grid_line_alpha = 0
p3.ygrid.band_fill_color = "olive"
p3.ygrid.band_fill_alpha = 0.1
p3.yaxis.axis_label = "Indivíduos"
p3.xaxis.axis_label = "Dias"

p4.grid.grid_line_alpha = 0
p4.ygrid.band_fill_color = "olive"
p4.ygrid.band_fill_alpha = 0.1
p4.yaxis.axis_label = "Indivíduos"
p4.xaxis.axis_label = "Dias"

# Incluindo as curvas
for data in saved_prediction[20:]:
    p3.line(pred_t, data[1],
           legend_label="Previsão Infectados",
           line_cap="round", line_dash = 'dashed',
           line_width=4, color="#42a5f5", line_alpha = 0.1)

    p3.line(pred_t, data[2],
           legend_label="Previsão Recuperados",
           line_cap="round", line_dash = 'dashed', line_width=4,
           color="#9c27b0", line_alpha = 0.07)


p3.line(td, Id,
       legend_label="Infectados",
       line_cap="round", line_width=3, color="#005cb2")

p3.line(td, Rd,
           legend_label="Recuperados",
           line_cap="round", line_width=3, color="#5e35b1")

if plot_all:
    for data in saved_prediction:
        p4.line(pred_t, data[0] + N*(1-saved_param['pop'][-1]),
               legend_label="Previsão Suscetiveis",
               line_cap="round", line_dash = 'dashed',
               line_width=4, color="#ff5722", line_alpha = 0.07)
    p4.line(td, Sd,
           legend_label="Suscetiveis",
           line_cap="round", line_width=3, color="#b71c1c")

# Colocando as legendas
p3.legend.click_policy="hide"
p3.legend.location = "top_left"

p4.legend.click_policy="hide"
p4.legend.location = "top_right"

show(column(p3,p4))

Análise de variação do \(R_0\)

[33]:
x = range(start_day, len(I))
beta_norm = [b*p for b,p in zip(saved_param['beta'],saved_param['pop'])]
Ro = [b*p/r for b,r,p in zip(saved_param['beta'],saved_param['r'],saved_param['pop'])]


fig1 = go.Figure()

fig1.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=beta_norm[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="β = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=beta_norm[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="β = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[beta_norm[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um β = %{y:.4f}."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[beta_norm[-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para β = %{y:.4f}."))
fig1.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.10, y=1.05),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro β',
                  title_text="Estimativa do parâmetro β - 🇨🇳 " + country['data']['name'])

fig1.show(renderer="png")
_images/notebooks_SIR_Real_Data_-_COVID_-_China_Model_Analyser_25_0.png
[32]:

fig2 = go.Figure()

fig2.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=saved_param['r'][:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig2.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=saved_param['r'][peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig2.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[saved_param['r'][peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um r = %{y:.4f}."))
fig2.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[saved_param['r'][-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para r = %{y:.4f}."))
fig2.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.07, y=1.06),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro r',
                  title_text="Estimativa do parâmetro r - 🇨🇳 " + country['data']['name'])

fig2.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_China_Model_Analyser_26_0.png
[31]:

fig3 = go.Figure()

fig3.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=Ro[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig3.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=Ro[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="Ro = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig3.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[Ro[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um Ro = %{y:.4f}."))
fig3.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[Ro[-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para Ro = %{y:.4f}."))
fig3.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.07, y=1.06),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro Ro',
                  title_text="Estimativa do parâmetro Ro - 🇨🇳 " + country['data']['name'])

fig3.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_China_Model_Analyser_27_0.png

Análise de confiança

[20]:

from PyAstronomy import pyasl

dI = np.gradient(pyasl.smooth(I, 13, "hamming"))
t = np.linspace(0, len(dI), len(dI))

signal = np.array([di >= 0 for di in dI[::-1]])


import plotly.graph_objects as go

fig = go.Figure()
fig.add_trace(go.Scatter(
                    x=t[signal],
                    y=dI[::-1][signal],
                    mode='lines',
                    name='Antes do pico - Derivada positiva',
                    line_shape='spline',
                    line = dict(color='#512da8', width=3)))
fig.add_trace(go.Scatter(
                    x=t[~signal],
                    y=dI[::-1][~signal],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico - Derivada negativa',
                    line = dict(color='#d8576b', width=3)))

fig.update_layout(template='xgridoff',
                  xaxis_title='Dias (invertido)',
                  yaxis_title='Taxa de variação de indivíduos',
                  title_text="Derivada da curva de Infectados - " + country['data']['name'])
# 'ggplot2', 'seaborn', 'simple_white', 'plotly',
# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',
# 'ygridoff', 'gridon', 'none'

fig.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_China_Model_Analyser_29_0.png
[22]:

peak_pos = np.argmax(~signal[::-1])
print("O pico da epidemia acontece no dia", peak_pos)

final_day = int(td[-1])

estimated_peaks = []
for data in saved_prediction:
    dI = np.gradient(data[1][:final_day])
    signal_pred = np.array([di >= 0 for di in dI[::-1]])
    estimated_peaks.append(len(Id) - np.argmax(signal_pred))
estimated_peaks = np.array(estimated_peaks)

O pico da epidemia acontece no dia 28
[23]:

import plotly.graph_objects as go

peak_error = np.abs(estimated_peaks - peak_pos)

fig1 = go.Figure()
fig1.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=peak_error[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="Erro de %{y} dias, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    x=td[peak_pos:],
                    y=peak_error[peak_pos-start_day:],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="Erro de %{y} dias, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[peak_error[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x} depois do começo da epidemia."))

fig1.update_layout(yaxis_range=[-1,31],
                  template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.20, y=1.0),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Erro (em dias) da estimativa do pico',
                  title_text="Erro da estimativa do pico - " + country['data']['name'])
# 'ggplot2', 'seaborn', 'simple_white', 'plotly',
# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',
# 'ygridoff', 'gridon', 'none'

fig1.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_China_Model_Analyser_31_0.png

SIR 🇮🇹 Itália - Dados Reais

Os dados

Para coletar os dados de outros países, foram utilizadas as APIs dos sites:

[1]:

import requests
import pandas as pd

covid_api = 'https://corona-api.com/countries/'
rest_countries = 'https://restcountries.eu/rest/v2/alpha/'
country = 'IT' # Alpha-2 ISO3166

data_json =  requests.get(covid_api + country).json()
country = requests.get(covid_api + country).json()

N = country['data']['population']

print(country['data']['name'])

Italy

Organizando os dados

[2]:

from datetime import datetime

df = pd.DataFrame(data_json['data']['timeline'])
df = df.sort_values('date').reset_index()

from datetime import datetime, timedelta
df['date'] = [datetime.fromisoformat(f) for f in df['date']]
df = df.drop_duplicates(subset='date', keep = 'last')

# Criando o vetor de tempo
first_date = df['date'].iloc[0]
size_days = (df['date'].iloc[-1] - df['date'].iloc[0]).days
date_vec = [first_date + timedelta(days=k) for k in range(size_days)]

new_df = pd.DataFrame(date_vec, columns=['date'])
new_df = pd.merge(new_df, df, how='left', on= 'date')
new_df = new_df.drop(columns= ['index',  'updated_at', 'is_in_progress'])

for col in new_df.columns[1:]:
    new_df[col] = new_df[col].interpolate(method='polynomial', order=1)
df = new_df.dropna()

df.head()

[2]:
date deaths confirmed active recovered new_confirmed new_recovered new_deaths
0 2020-01-30 0.0 2.000000 2.000000 0.0 2.000000 0.0 0.0
1 2020-01-31 0.0 2.000000 2.000000 0.0 0.000000 0.0 0.0
2 2020-02-01 0.0 2.142857 2.142857 0.0 0.142857 0.0 0.0
3 2020-02-02 0.0 2.285714 2.285714 0.0 0.285714 0.0 0.0
4 2020-02-03 0.0 2.428571 2.428571 0.0 0.428571 0.0 0.0

Visualizando os dados

[3]:

from bokeh.models   import Legend, ColumnDataSource, RangeTool, LinearAxis, Range1d, HoverTool
from bokeh.palettes import brewer, Inferno256
from bokeh.plotting import figure, show
from bokeh.layouts  import column
from bokeh.io       import output_notebook

output_notebook()

import numpy as np

# Criando os valores para legenda no plot
year =  [str(int(d.year)) for d in df['date'] ]
month = [("0"+str(int(d.month)))[-2:] for d in df['date'] ]
day =   [("0"+str(int(d.day)))[-2:] for d in df['date'] ]

# Criando a fonte de dados
source = ColumnDataSource(data={
    'Data'       : df['date'].values,
    'd': day, 'm': month, 'y': year,
    'Infectados Acc' : df['confirmed'].values,
    'Mortes'     : df['deaths'].values,
    'Ativo'     : df['active'].values,
    'Recuperados': df['recovered'].values
})


# Criando a figura
p = figure(plot_height=500,
           plot_width=600,
           x_axis_type="datetime",
           tools="",
           #y_axis_type="log",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

# Incluindo as curvas
i_p = p.line(x='Data', y='Ativo',
             legend_label="Infectados",
             line_cap="round", line_width=5, color="#c62828", source=source)
m_p = p.line(x='Data', y='Mortes',
             legend_label="Mortes",
             line_cap="round", line_width=5, color="#512da8", source=source)
c_p = p.line(x='Data', y='Infectados Acc',
             legend_label="Infectados Acc",
             line_cap="round", line_width=5, color="#0288d1", source=source)
r_p = p.line(x='Data', y='Recuperados',
             legend_label="Recuperados",
             line_cap="round", line_width=5, color="#388e3c", source=source)

# Colocando as legendas
p.legend.click_policy="hide"
# p.legend.location = "top_left"
p.legend.location = "top_left"

# Incluindo a ferramenta de hover
p.add_tools(HoverTool(
    tooltips=[
        ( 'Indivíduos', '$y{i}'),
        ( 'Data',       '@d/@m/@y' ),
    ],
    renderers=[
        m_p, i_p, c_p, r_p
    ]
))

show(p)

Loading BokehJS ...
Verificando os dados
[4]:

dif_I = np.diff(df['active'])
cum = []
cum.append(dif_I[0])

for k, i in enumerate(dif_I):
    cum.append(cum[-1] + i)

cum = np.array(cum)
cum += df['deaths'].to_numpy() + df['recovered'].to_numpy()

print("Erro entre casos acumulados e valores de confirmados: {}".format(
        round(sum((cum - df['confirmed'].values)**2 / len(cum)),2) ) )

Erro entre casos acumulados e valores de confirmados: 4.0
Criando os dados SIR
[5]:

I = df['active'].to_numpy()
R = df['recovered'].to_numpy()
M = df['deaths'].to_numpy()
S = N - R - I

# Creating the time vector
t = np.linspace(0, len(I), len(I))


Sd, Id, Md, Rd, td = S, I, M, R, t

Estimando utilizando todos os dados

[6]:

from models import *

dataset = dict(S=Sd, I=Id, R=Rd)

# Create the model
sir_model = ss.SIR(pop=N, focus=["S", "I", "R"])

# Adjust the parameters
sir_model.fit(dataset, td,
              search_pop=True,
              pop_sens=[0.001,0.01],
              beta_sens=[100,100],
              r_sens=[10000,100])

# Predict the model
sim_res = sir_model.predict((Sd[0],Id[0], Rd[0]), td)

Loading BokehJS ...
         ├─ S(0) ─ I(0) ─ R(0) ─  [60340326.0, 2.0, 0.0]
         ├─ beta ─   1   r ─   0.14285714285714285
         ├─ beta bound ─   0.01  ─  100
         ├─ r bound ─   1.4285714285714285e-05  ─  14.285714285714285
         ├─ equation weights ─   [1.6598704101560237e-08, 1.9910826361712874e-05, 2.2503737227789617e-05]
         ├─ Running on ─  differential_evolution SciPy Search Algorithm
         └─ Defined at:  65.58644881554696  ─  0.023308094571579452

[7]:

print("Parâmetros estimados: ", sir_model.parameters)
print("Suposto Ro: ", sir_model.parameters[0] * sir_model.parameters[-1] / sir_model.parameters[1])
print("Dias contaminados: ", 1 / sir_model.parameters[1])

Parâmetros estimados:  [6.55864488e+01 2.33080946e-02 3.10705183e-03]
Suposto Ro:  8.742906691463745
Dias contaminados:  42.90354996325364
[56]:

p = figure(plot_height=500,
           plot_width=700,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

p.line(t, I,
       legend_label="Infectados", color="#ff6659", line_width=4)
p.line(t, R,
       legend_label="Removidos", color="#76d275", line_width=4)

# Show the results
p.line(td, sim_res[1],
       legend_label="Infectados - Modelo", line_dash="dashed", color="#d32f2f", line_width=3)
p.line(td, sim_res[2],
    legend_label="Removidos - Modelo", line_dash="dashed", color="#43a047", line_width=3)
p.line(td, sim_res[0],
      legend_label="Suscetíveis Ponderados - Modelo", line_dash="dashed", color="#1e88e5", line_width=3)
p.add_layout(p.legend[0], 'right')

show(p)

Monte Carlo

Nesta parte, faremos um teste aumentando a quantidade de amostras de treinamento e prevendo o momento do pico da epidemia a medida que mais dias são utilizados para treinamento. Esse estudo vai possibilitar a análise da certeza da previsão do pico da epidemia antes desse acontecer.

[13]:

saved_param = {'r':[], 'beta':[], 'pop':[]}
saved_prediction = []

start_day = 35
pred_t = np.array(range(int(td[-1]) + 120))

for i in range(start_day,len(I),1):

    dataset = dict(S=Sd[:i], I=Id[:i], R=Rd[:i])
    td_ = td[:i]

    # Create the model
    sir_model = ss.SIR(pop=N, focus=["S", "I", "R"], verbose = False)

    # Adjust the parameters
    sir_model.fit(dataset, td_,
                  search_pop=True,
                  pop_sens=[0.001,0.01],
                  beta_sens=[100000,100],
                  r_sens=[100000,100])

    saved_param['beta'].append(sir_model.parameters[0])
    saved_param['r'].append(sir_model.parameters[1])
    saved_param['pop'].append(sir_model.parameters[2])

    saved_prediction.append(sir_model.predict((Sd[0],Id[0],Rd[0]), pred_t))
[27]:
import pickle

with open("./Italy_mc_runs.pickle", "wb") as handle:
    pickle.dump({"pars":saved_param, "pred":saved_prediction}, handle)

Análise do uso do sistema de saúde

Nesta análise, mostramos o erro percentual do quanto antes do pico, conseguimos prever a quantidade de pessoas que realmente serão identificadas como infectadas, uma vez que \(R(\infty)\) é a quantidade de pessoas recuperadas totais, daquelas que foram noficadas como infectadas no sistema de saúde. Desta forma segue o erro proporcional do erro a medida em que novos dados diários foram incluidos no modelo:

[71]:

x = range(start_day, len(I))

usage_error = [ 100 * abs(p*N - Rd[-1]) / Rd[-1] for p in saved_param['pop']]

fig9 = go.Figure()

fig9.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=usage_error[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="ε(%) = %{y:.0f}, <br> com %{x:.0f} dias de dados."))
fig9.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=usage_error[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="ε(%) = %{y:.0f}, <br> com %{x:.0f} dias de dados."))
fig9.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[usage_error[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um ε(%) = %{y:.0f}."))
fig9.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[usage_error[-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para ε(%) = %{y:.4f}."))
fig9.update_layout(template='xgridoff', yaxis_range=[-1,100],
                  legend_orientation="h", legend=dict(x=0.10, y=1.05),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Erro (%) do R(∞) - notificados',
                  title_text="Erro da previsão do uso do sistema de saúde - 🇮🇹 " + country['data']['name'])

fig9.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Italy_Model_Analyser_18_0.png

Visualizando as previsões de \(I(t)\)

Vamos analisar as previsões quando somente os dados antes do pico são fornecidos ao modelo, e as previões utilizando os dados após o pico:

[53]:

p2 = figure(plot_height=500,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p2.grid.grid_line_alpha = 0
p2.ygrid.band_fill_color = "olive"
p2.ygrid.band_fill_alpha = 0.1
p2.yaxis.axis_label = "Indivíduos"
p2.xaxis.axis_label = "Dias"

# Incluindo as curvas
p2.line(td, Id,
       legend_label="Infectados",
       line_cap="round", line_width=3, color="#c62828")

for data in saved_prediction[45:]:
    p2.line(pred_t[:len(td)], -data[1][:len(td)],
           legend_label="Previsão Infectados - Depois do pico",
           line_cap="round", line_dash="dashed", line_width=4, color="#ffa000", line_alpha = 0.1)

for data in saved_prediction[:45]:
    p2.line(pred_t[:len(td)], -data[1][:len(td)],
           legend_label="Previsão Infectados - Antes do pico",
           line_cap="round", line_dash="dashed", line_width=4, color="#42a5f5", line_alpha = 0.1)

# Colocando as legendas
p2.legend.click_policy="hide"
p2.legend.location = "bottom_left"

show(p2)

Visualizando os ajustes dos grupos

Aqui vamos analisar as previsões obtidas para cada grupo \(S(t)\), \(I(t)\) e \(R(t)\), a medida que mais dias foram informados ao modelo.

[54]:

p3 = figure(plot_height=350,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - I(t) e R(t) - " + country['data']['name'])

p4 = figure(plot_height=350,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - S(t) - " + country['data']['name'])

plot_all = True

# Preparando o estilo
p3.grid.grid_line_alpha = 0
p3.ygrid.band_fill_color = "olive"
p3.ygrid.band_fill_alpha = 0.1
p3.yaxis.axis_label = "Indivíduos"
p3.xaxis.axis_label = "Dias"

p4.grid.grid_line_alpha = 0
p4.ygrid.band_fill_color = "olive"
p4.ygrid.band_fill_alpha = 0.1
p4.yaxis.axis_label = "Indivíduos"
p4.xaxis.axis_label = "Dias"

# Incluindo as curvas
for data in saved_prediction[10:]:
    p3.line(pred_t, data[1],
           legend_label="Previsão Infectados",
           line_cap="round", line_dash = 'dashed',
           line_width=4, color="#42a5f5", line_alpha = 0.1)

    p3.line(pred_t, data[2],
           legend_label="Previsão Recuperados",
           line_cap="round", line_dash = 'dashed', line_width=4,
           color="#9c27b0", line_alpha = 0.07)


p3.line(td, Id,
       legend_label="Infectados",
       line_cap="round", line_width=3, color="#005cb2")

p3.line(td, Rd,
           legend_label="Recuperados",
           line_cap="round", line_width=3, color="#5e35b1")

if plot_all:
    for data in saved_prediction[10:]:
        p4.line(pred_t, data[0] + N*(1-saved_param['pop'][-1]),
               legend_label="Previsão Suscetiveis",
               line_cap="round", line_dash = 'dashed',
               line_width=4, color="#ff5722", line_alpha = 0.07)
    p4.line(td, Sd,
           legend_label="Suscetiveis",
           line_cap="round", line_width=3, color="#b71c1c")

# Colocando as legendas
p3.legend.click_policy="hide"
p3.legend.location = "top_left"

p4.legend.click_policy="hide"
p4.legend.location = "top_right"

show(column(p3,p4))

Análise de variação do \(R_0\)

[67]:


x = range(start_day, len(I))
beta_norm = [b*p for b,p in zip(saved_param['beta'],saved_param['pop'])]
Ro = [b*p/r for b,r,p in zip(saved_param['beta'],saved_param['r'],saved_param['pop'])]


fig1 = go.Figure()

fig1.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=beta_norm[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="β = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=beta_norm[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="β = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[beta_norm[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um β = %{y:.4f}."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[beta_norm[-1]],
                    marker_symbol="triangle-right-dot", name="Valor de convergência",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> converge para β = %{y:.4f}."))
fig1.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.10, y=1.05),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro β',
                  title_text="Estimativa do parâmetro β - 🇮🇹 " + country['data']['name'])

fig1.show(renderer="png")
_images/notebooks_SIR_Real_Data_-_COVID_-_Italy_Model_Analyser_24_0.png
[68]:

fig2 = go.Figure()

fig2.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=saved_param['r'][:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig2.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=saved_param['r'][peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig2.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[saved_param['r'][peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um r = %{y:.4f}."))
fig2.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[saved_param['r'][-1]],
                    marker_symbol="triangle-right-dot", name="Valor de convergência",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> converge para r = %{y:.4f}."))
fig2.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.07, y=1.06),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro r',
                  title_text="Estimativa do parâmetro r - 🇮🇹 " + country['data']['name'])

fig2.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Italy_Model_Analyser_25_0.png
[69]:

fig3 = go.Figure()

fig3.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=Ro[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig3.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=Ro[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="Ro = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig3.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[Ro[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um Ro = %{y:.4f}."))
fig3.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[Ro[-1]],
                    marker_symbol="triangle-right-dot", name="Valor de atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para Ro = %{y:.4f}."))
fig3.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.07, y=1.06),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro Ro',
                  title_text="Estimativa do parâmetro Ro - 🇮🇹 " + country['data']['name'])

fig3.show(renderer="png")


_images/notebooks_SIR_Real_Data_-_COVID_-_Italy_Model_Analyser_26_0.png

Análise de confiança

[70]:

from PyAstronomy import pyasl

dI = np.gradient(pyasl.smooth(I, 13, "hamming"))
t = np.linspace(0, len(dI), len(dI))

signal = np.array([di >= 0 for di in dI[::-1]])


import plotly.graph_objects as go

fig = go.Figure()
fig.add_trace(go.Scatter(
                    x=t[signal],
                    y=dI[::-1][signal],
                    mode='lines',
                    name='Antes do pico - Derivada positiva',
                    line_shape='spline',
                    line = dict(color='#512da8', width=3)))
fig.add_trace(go.Scatter(
                    x=t[~signal],
                    y=dI[::-1][~signal],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico - Derivada negativa',
                    line = dict(color='#d8576b', width=3)))

fig.update_layout(template='xgridoff',
                  xaxis_title='Dias (invertido)',
                  yaxis_title='Taxa de variação de indivíduos',
                  title_text="Derivada da curva de Infectados - 🇮🇹 " + country['data']['name'])
# 'ggplot2', 'seaborn', 'simple_white', 'plotly',
# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',
# 'ygridoff', 'gridon', 'none'

fig.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Italy_Model_Analyser_28_0.png
[22]:

peak_pos = len(Id) - np.argmax(signal)
print("O pico da epidemia acontece no dia", peak_pos)

start_day = 35
final_day = int(td[-1])

estimated_peaks = []
for data in saved_prediction:
    dI = np.gradient(data[1][:final_day])
    signal_pred = np.array([di >= 0 for di in dI[::-1]])
    estimated_peaks.append(len(Id) - np.argmax(signal_pred))
estimated_peaks = np.array(estimated_peaks)

O pico da epidemia acontece no dia 82
[47]:

import plotly.graph_objects as go

peak_error = np.abs(estimated_peaks - peak_pos)

fig1 = go.Figure()
fig1.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=peak_error[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="Erro de %{y} dias, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    x=td[peak_pos:],
                    y=peak_error[peak_pos-start_day:],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="Erro de %{y} dias, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[peak_error[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x} depois do começo da epidemia."))

fig1.update_layout(yaxis_range=[-1,31], template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.20, y=1.0),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Erro (em dias) da estimativa do pico',
                  title_text="Erro da estimativa do pico - 🇮🇹 " + country['data']['name'])
# 'ggplot2', 'seaborn', 'simple_white', 'plotly',
# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',
# 'ygridoff', 'gridon', 'none'

fig1.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Italy_Model_Analyser_30_0.png

SIR 🇩🇪 Alemanha - Dados Reais

Os dados

Para coletar os dados de outros países, foram utilizadas as APIs dos sites:

[1]:

import requests
import pandas as pd

covid_api = 'https://corona-api.com/countries/'
rest_countries = 'https://restcountries.eu/rest/v2/alpha/'
country = 'DE' # Alpha-2 ISO3166

data_json =  requests.get(covid_api + country).json()
country = requests.get(covid_api + country).json()

N = country['data']['population']

print(country['data']['name'])

Germany

Organizando os dados

[2]:

from datetime import datetime

df = pd.DataFrame(data_json['data']['timeline'])
df = df.sort_values('date').reset_index()

from datetime import datetime, timedelta
df['date'] = [datetime.fromisoformat(f) for f in df['date']]
df = df.drop_duplicates(subset='date', keep = 'last')

# Criando o vetor de tempo
first_date = df['date'].iloc[0]
size_days = (df['date'].iloc[-1] - df['date'].iloc[0]).days
date_vec = [first_date + timedelta(days=k) for k in range(size_days)]

new_df = pd.DataFrame(date_vec, columns=['date'])
new_df = pd.merge(new_df, df, how='left', on= 'date')
new_df = new_df.drop(columns= ['index',  'updated_at', 'is_in_progress'])

for col in new_df.columns[1:]:
    new_df[col] = new_df[col].interpolate(method='polynomial', order=1)
df = new_df.dropna()

df.head()

[2]:
date deaths confirmed active recovered new_confirmed new_recovered new_deaths
0 2020-01-27 0.0 4.0 4.0 0.0 4.0 0.0 0.0
1 2020-01-28 0.0 4.0 4.0 0.0 0.0 0.0 0.0
2 2020-01-29 0.0 4.0 4.0 0.0 0.0 0.0 0.0
3 2020-01-30 0.0 5.0 5.0 0.0 1.0 0.0 0.0
4 2020-01-31 0.0 8.0 8.0 0.0 3.0 0.0 0.0

Visualizando os dados

[3]:

from bokeh.models   import Legend, ColumnDataSource, RangeTool, LinearAxis, Range1d, HoverTool
from bokeh.palettes import brewer, Inferno256
from bokeh.plotting import figure, show
from bokeh.layouts  import column
from bokeh.io       import output_notebook

output_notebook()

import numpy as np

# Criando os valores para legenda no plot
year =  [str(int(d.year)) for d in df['date'] ]
month = [("0"+str(int(d.month)))[-2:] for d in df['date'] ]
day =   [("0"+str(int(d.day)))[-2:] for d in df['date'] ]

# Criando a fonte de dados
source = ColumnDataSource(data={
    'Data'       : df['date'].values,
    'd': day, 'm': month, 'y': year,
    'Infectados Acc' : df['confirmed'].values,
    'Mortes'     : df['deaths'].values,
    'Ativo'     : df['active'].values,
    'Recuperados': df['recovered'].values
})


# Criando a figura
p = figure(plot_height=500,
           plot_width=600,
           x_axis_type="datetime",
           tools="",
           #y_axis_type="log",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

# Incluindo as curvas
i_p = p.line(x='Data', y='Ativo',
             legend_label="Infectados",
             line_cap="round", line_width=5, color="#c62828", source=source)
m_p = p.line(x='Data', y='Mortes',
             legend_label="Mortes",
             line_cap="round", line_width=5, color="#512da8", source=source)
c_p = p.line(x='Data', y='Infectados Acc',
             legend_label="Infectados Acc",
             line_cap="round", line_width=5, color="#0288d1", source=source)
r_p = p.line(x='Data', y='Recuperados',
             legend_label="Recuperados",
             line_cap="round", line_width=5, color="#388e3c", source=source)

# Colocando as legendas
p.legend.click_policy="hide"
# p.legend.location = "top_left"
p.legend.location = "top_left"

# Incluindo a ferramenta de hover
p.add_tools(HoverTool(
    tooltips=[
        ( 'Indivíduos', '$y{i}'),
        ( 'Data',       '@d/@m/@y' ),
    ],
    renderers=[
        m_p, i_p, c_p, r_p
    ]
))

show(p)

Loading BokehJS ...
Verificando os dados
[4]:

dif_I = np.diff(df['active'])
cum = []
cum.append(dif_I[0])

for k, i in enumerate(dif_I):
    cum.append(cum[-1] + i)

cum = np.array(cum)
cum += df['deaths'].to_numpy() + df['recovered'].to_numpy()

print("Erro entre casos acumulados e valores de confirmados: {}".format(
        round(sum((cum - df['confirmed'].values)**2 / len(cum)),2) ) )

Erro entre casos acumulados e valores de confirmados: 16.0
Criando os dados SIR
[5]:

I = df['active'].to_numpy()
R = df['recovered'].to_numpy()
M = df['deaths'].to_numpy()
S = N - R - I

# Creating the time vector
t = np.linspace(0, len(I), len(I))


Sd, Id, Md, Rd, td = S, I, M, R, t

Estimando utilizando todos os dados

[6]:

from models import *

dataset = dict(S=Sd, I=Id, R=Rd)

# Create the model
sir_model = ss.SIR(pop=N, focus=["S", "I", "R"])

# Adjust the parameters
sir_model.fit(dataset, td,
              search_pop=True,
              pop_sens=[0.001,0.05],
              beta_sens=[100,1000],
              r_sens=[100,100])


# Predict the model
sim_res = sir_model.predict((Sd[0],Id[0],Rd[0]), td)

Loading BokehJS ...
         ├─ S(0) ─ I(0) ─ R(0) ─  [81802253.0, 4.0, 0.0]
         ├─ beta ─   1   r ─   0.14285714285714285
         ├─ beta bound ─   0.01  ─  1000
         ├─ r bound ─   0.0014285714285714286  ─  14.285714285714285
         ├─ equation weights ─   [1.223645216108376e-08, 4.8013887093497817e-05, 1.7124646680667163e-05]
         ├─ Running on ─  differential_evolution SciPy Search Algorithm
         └─ Defined at:  105.13818075597061  ─  0.05759270422192064

[7]:

print("Parâmetros estimados: ", sir_model.parameters)
print("Suposto Ro: ", sir_model.parameters[0] * sir_model.parameters[-1] / sir_model.parameters[1])
print("Dias contaminados: ", 1 / sir_model.parameters[1])

Parâmetros estimados:  [1.05138181e+02 5.75927042e-02 2.08772990e-03]
Suposto Ro:  3.811248776553262
Dias contaminados:  17.36331039686421
[19]:

p = figure(plot_height=500,
           plot_width=800,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"

p.line(t, I,
       legend_label="Infectados", color="#ff6659", line_width=4)
p.line(t, R,
       legend_label="Removidos", color="#76d275", line_width=4)

# Show the results
p.line(td, sim_res[1],
       legend_label="Infectados - Modelo", line_dash="dashed", color="#d32f2f", line_width=3)
p.line(td, sim_res[2],
    legend_label="Removidos - Modelo", line_dash="dashed", color="#43a047", line_width=3)
p.line(td, sim_res[0],
      legend_label="Suscetíveis Ponderados - Modelo", line_dash="dashed", color="#1e88e5", line_width=3)
p.add_layout(p.legend[0], 'right')

show(p)

Monte Carlo

Nesta parte, faremos um teste aumentando a quantidade de amostras de treinamento e prevendo o momento do pico da epidemia a medida que mais dias são utilizados para treinamento. Esse estudo vai possibilitar a análise da certeza da previsão do pico da epidemia antes desse acontecer.

[12]:

saved_param = {'r':[], 'beta':[], 'pop':[]}
saved_prediction = []

start_day = 45
pred_t = np.array(range(int(td[-1]) + 120))

for i in range(start_day, len(I)):

    dataset = dict(S=Sd[:i], I=Id[:i], R=Rd[:i])
    td_ = td[:i]

    # Create the model
    sir_model = ss.SIR(pop=N, focus=["S", "I", "R"], verbose = False)

    # Adjust the parameters
    sir_model.fit(dataset, td_,
                  search_pop=True,
                  pop_sens=[0.001,0.04],
                  beta_sens=[100,1000],
                  r_sens=[100,100])

    saved_param['beta'].append(sir_model.parameters[0])
    saved_param['r'].append(sir_model.parameters[1])
    saved_param['pop'].append(sir_model.parameters[2])

    saved_prediction.append(sir_model.predict((Sd[0],Id[0], Rd[0]), pred_t))
[13]:
import pickle

with open("./Germany_mc_runs.pickle", "wb") as handle:
    pickle.dump({"pars":saved_param, "pred":saved_prediction}, handle)

Análise do uso do sistema de saúde

Nesta análise, mostramos o erro percentual do quanto antes do pico, conseguimos prever a quantidade de pessoas que realmente serão identificadas como infectadas, uma vez que \(R(\infty)\) é a quantidade de pessoas recuperadas totais, daquelas que foram noficadas como infectadas no sistema de saúde. Desta forma segue o erro proporcional do erro a medida em que novos dados diários foram incluidos no modelo:

[41]:

x = range(start_day, len(I))
usage_error = [ 100 * abs(p*N - Rd[-1]) / Rd[-1] for p in saved_param['pop']]

fig9 = go.Figure()

fig9.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=usage_error[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="ε(%) = %{y:.0f}, <br> com %{x:.0f} dias de dados."))
fig9.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=usage_error[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="ε(%) = %{y:.0f}, <br> com %{x:.0f} dias de dados."))
fig9.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[usage_error[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um ε(%) = %{y:.0f}."))
fig9.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[usage_error[-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para ε(%) = %{y:.4f}."))
fig9.update_layout(template='xgridoff', yaxis_range=[-1,100],
                  legend_orientation="h", legend=dict(x=0.10, y=1.05),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Erro (%) do R(∞) - notificados',
                  title_text="Erro da previsão do uso do sistema de saúde - 🇩🇪 " + country['data']['name'])

fig9.show(renderer="png")


_images/notebooks_SIR_Real_Data_-_COVID_-_Germany_Model_Analyser_18_0.png

Visualizando as previsões de \(I(t)\)

Vamos analisar as previsões quando somente os dados antes do pico são fornecidos ao modelo, e as previões utilizando os dados após o pico:

[17]:

visual_peak = 70

p2 = figure(plot_height=500,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - " + country['data']['name'])

# Preparando o estilo
p2.grid.grid_line_alpha = 0
p2.ygrid.band_fill_color = "olive"
p2.ygrid.band_fill_alpha = 0.1
p2.yaxis.axis_label = "Indivíduos"
p2.xaxis.axis_label = "Dias"

# Incluindo as curvas
p2.line(td, Id,
       legend_label="Infectados",
       line_cap="round", line_width=3, color="#c62828")

for data in saved_prediction[visual_peak-start_day:]:
    p2.line(pred_t[:len(td)], -data[1][:len(td)],
           legend_label="Previsão Infectados - Depois do pico",
           line_cap="round", line_dash="dashed", line_width=4, color="#ffa000", line_alpha = 0.1)

for data in saved_prediction[20:visual_peak-start_day]:
    p2.line(pred_t[:len(td)], -data[1][:len(td)],
           legend_label="Previsão Infectados - Antes do pico",
           line_cap="round", line_dash="dashed", line_width=4, color="#42a5f5", line_alpha = 0.1)

# Colocando as legendas
p2.legend.click_policy="hide"
p2.legend.location = "bottom_left"

show(p2)

Visualizando os ajustes dos grupos

Aqui vamos analisar as previsões obtidas para cada grupo \(S(t)\), \(I(t)\) e \(R(t)\), a medida que mais dias foram informados ao modelo.

[33]:

p3 = figure(plot_height=350,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - I(t) e R(t) - " + country['data']['name'])

p4 = figure(plot_height=350,
           plot_width=600,
           tools="",
           toolbar_location=None,
           title="Evolução do COVID - S(t) - " + country['data']['name'])

plot_all = True

# Preparando o estilo
p3.grid.grid_line_alpha = 0
p3.ygrid.band_fill_color = "olive"
p3.ygrid.band_fill_alpha = 0.1
p3.yaxis.axis_label = "Indivíduos"
p3.xaxis.axis_label = "Dias"

p4.grid.grid_line_alpha = 0
p4.ygrid.band_fill_color = "olive"
p4.ygrid.band_fill_alpha = 0.1
p4.yaxis.axis_label = "Indivíduos"
p4.xaxis.axis_label = "Dias"

# Incluindo as curvas
for data in saved_prediction[20:]:
    p3.line(pred_t, data[1],
           legend_label="Previsão Infectados",
           line_cap="round", line_dash = 'dashed',
           line_width=4, color="#42a5f5", line_alpha = 0.1)

    p3.line(pred_t, data[2],
           legend_label="Previsão Recuperados",
           line_cap="round", line_dash = 'dashed', line_width=4,
           color="#9c27b0", line_alpha = 0.07)


p3.line(td, Id,
       legend_label="Infectados",
       line_cap="round", line_width=3, color="#005cb2")

p3.line(td, Rd,
           legend_label="Recuperados",
           line_cap="round", line_width=3, color="#5e35b1")

if plot_all:
    for data in saved_prediction:
        p4.line(pred_t, data[0] + N*(1-saved_param['pop'][-1]),
               legend_label="Previsão Suscetiveis",
               line_cap="round", line_dash = 'dashed',
               line_width=4, color="#ff5722", line_alpha = 0.07)
    p4.line(td, Sd,
           legend_label="Suscetiveis",
           line_cap="round", line_width=3, color="#b71c1c")

# Colocando as legendas
p3.legend.click_policy="hide"
p3.legend.location = "top_left"

p4.legend.click_policy="hide"
p4.legend.location = "top_right"

show(column(p3,p4))

Análise de variação do \(R_0\)

[42]:

x = range(start_day, len(I))
beta_norm = [b*p for b,p in zip(saved_param['beta'],saved_param['pop'])]
Ro = [b*p/r for b,r,p in zip(saved_param['beta'],saved_param['r'],saved_param['pop'])]


fig1 = go.Figure()

fig1.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=beta_norm[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="β = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=beta_norm[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="β = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[beta_norm[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um β = %{y:.4f}."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[beta_norm[-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para β = %{y:.4f}."))
fig1.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.10, y=1.05),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro β',
                  title_text="Estimativa do parâmetro β - 🇩🇪 " + country['data']['name'])

fig1.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Germany_Model_Analyser_24_0.png
[43]:

fig2 = go.Figure()

fig2.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=saved_param['r'][:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig2.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=saved_param['r'][peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig2.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[saved_param['r'][peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um r = %{y:.4f}."))
fig2.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[saved_param['r'][-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para r = %{y:.4f}."))
fig2.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.07, y=1.06),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro r',
                  title_text="Estimativa do parâmetro r - 🇩🇪 " + country['data']['name'])

fig2.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Germany_Model_Analyser_25_0.png
[44]:

fig3 = go.Figure()

fig3.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=Ro[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="r = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig3.add_trace(go.Scatter(
                    x=td[peak_pos:-1],
                    y=Ro[peak_pos-start_day:-1],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="Ro = %{y:.4f}, <br> com %{x:.0f} dias de dados."))
fig3.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[Ro[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x}, com um Ro = %{y:.4f}."))
fig3.add_trace(go.Scatter(
                    mode="markers", x=[td[-1]], y=[Ro[-1]],
                    marker_symbol="triangle-right-dot", name="Valor atual",
                    marker_line_color="#a00037", marker_color="#ff5c8d",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="No dia %{x} da epidemia, <br> convergindo para Ro = %{y:.4f}."))
fig3.update_layout(template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.07, y=1.06),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Parâmetro Ro',
                  title_text="Estimativa do parâmetro Ro - 🇩🇪 " + country['data']['name'])

fig3.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Germany_Model_Analyser_26_0.png

Análise de confiança

[22]:

from PyAstronomy import pyasl

dI = np.gradient(pyasl.smooth(I, 13, "hamming"))
t = np.linspace(0, len(dI), len(dI))

signal = np.array([di >= 0 for di in dI[::-1]])


import plotly.graph_objects as go

fig = go.Figure()
fig.add_trace(go.Scatter(
                    x=t[signal],
                    y=dI[::-1][signal],
                    mode='lines',
                    name='Antes do pico - Derivada positiva',
                    line_shape='spline',
                    line = dict(color='#512da8', width=3)))
fig.add_trace(go.Scatter(
                    x=t[~signal],
                    y=dI[::-1][~signal],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico - Derivada negativa',
                    line = dict(color='#d8576b', width=3)))

fig.update_layout(template='xgridoff',
                  xaxis_title='Dias (invertido)',
                  yaxis_title='Taxa de variação de indivíduos',
                  title_text="Derivada da curva de Infectados - " + country['data']['name'])
# 'ggplot2', 'seaborn', 'simple_white', 'plotly',
# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',
# 'ygridoff', 'gridon', 'none'

fig.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Germany_Model_Analyser_28_0.png
[23]:

peak_pos = len(Id) - np.argmax(signal)
print("O pico da epidemia acontece no dia", peak_pos)

final_day = int(td[-1])

estimated_peaks = []
for data in saved_prediction:
    dI = np.gradient(data[1][:final_day])
    signal_pred = np.array([di >= 0 for di in dI[::-1]])
    estimated_peaks.append(len(Id) - np.argmax(signal_pred))
estimated_peaks = np.array(estimated_peaks)

O pico da epidemia acontece no dia 71
[24]:

import plotly.graph_objects as go

peak_error = np.abs(estimated_peaks - peak_pos)

fig1 = go.Figure()
fig1.add_trace(go.Scatter(
                    x=td[start_day:peak_pos],
                    y=peak_error[:peak_pos-start_day],
                    mode='lines',
                    name='Antes do pico',
                    line_shape='spline',
                    line = dict(color='royalblue', width=3),
                    hovertemplate="Erro de %{y} dias, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    x=td[peak_pos:],
                    y=peak_error[peak_pos-start_day:],
                    mode='lines',
                    line_shape='spline',
                    name='Depois do pico',
                    line = dict(color='royalblue', width=3, dash='dot'),
                    hovertemplate="Erro de %{y} dias, <br> com %{x:.0f} dias de dados."))
fig1.add_trace(go.Scatter(
                    mode="markers", x=[peak_pos], y=[peak_error[peak_pos-start_day]],
                    marker_symbol="hexagram-dot", name="Momento do pico",
                    marker_line_color="midnightblue", marker_color="lightskyblue",
                    marker_line_width=2, marker_size=15,
                    hovertemplate="Pico no dia %{x} depois do começo da epidemia."))

fig1.update_layout( template='xgridoff',
                  legend_orientation="h", legend=dict(x=0.20, y=1.0),
                  xaxis_title='Dias (desde o começo da epidemia) utilizados para treinamento',
                  yaxis_title='Erro (em dias) da estimativa do pico',
                  title_text="Erro da estimativa do pico - " + country['data']['name'])
# 'ggplot2', 'seaborn', 'simple_white', 'plotly',
# 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff',
# 'ygridoff', 'gridon', 'none'

fig1.show(renderer="png")

_images/notebooks_SIR_Real_Data_-_COVID_-_Germany_Model_Analyser_30_0.png

epidemicModels

models package

Cost Functions

models.cost_functions.PrintException(e)
models.cost_functions.cost_NSIR(self, pars, dataset, initial, t, w)
models.cost_functions.cost_SEIR(self, pars, dataset, initial, t, w)

The function to compute the error to guide the learning algorithm. It computes the quadratic error.

Parameters:
  • pars (tuple) – Tuple with Beta and r parameters, respectivelly.
  • dataset (list) – The dataset with the respective S, I and R arrays.
  • initial (array) – The initial values of suceptible and infected, respectivelly.
  • t (array) – The time respective to each sample.
  • w (array) – The weight respective to the suceptible and infected errors.
Returns:

The sum of the quadratic error, between simulated and real data.

Return type:

float

models.cost_functions.cost_SIR(self, pars, dataset, initial, t, w)

The function to compute the error to guide the learning algorithm. It computes the quadratic error.

Parameters:
  • p (tuple) – Tuple with Beta and r parameters, respectivelly.
  • S (array) – The suceptible data values.
  • I (array) – The infected data values.
  • initial (array) – The initial values of suceptible and infected, respectivelly.
  • t (array) – The time respective to each sample.
  • w (array) – The weight respective to the suceptible and infected errors.
Returns:

The sum of the quadratic error, between simulated and real data.

Return type:

float

models.cost_functions.cost_SIRD(self, pars, dataset, initial, t, w)

The function to compute the error to guide the learning algorithm. It computes the quadratic error.

Parameters:
  • p (tuple) – Tuple with Beta and r parameters, respectivelly.
  • S (array) – The suceptible data values.
  • I (array) – The infected data values.
  • initial (array) – The initial values of suceptible and infected, respectivelly.
  • t (array) – The time respective to each sample.
  • w (array) – The weight respective to the suceptible and infected errors.
Returns:

The sum of the quadratic error, between simulated and real data.

Return type:

float

models.cost_functions.cost_dSIR(self, pars, dataset, initial, t, w)

The function to compute the error to guide the learning algorithm. It computes the quadratic error.

Parameters:
  • p (tuple) – Tuple with Beta and r parameters, respectivelly.
  • S (array) – The suceptible data values.
  • I (array) – The infected data values.
  • initial (array) – The initial values of suceptible and infected, respectivelly.
  • t (array) – The time respective to each sample.
  • w (array) – The weight respective to the suceptible and infected errors.
Returns:

The sum of the quadratic error, between simulated and real data.

Return type:

float

Differential Models

models.differential_models.NSIR(self, y, t, beta, r, betan, alpha, rn, *args)
models.differential_models.SEIR(self, y, t, Beta, r, sigma)

The function that computes the diferential set of equations of the SEIR Epidemic Model.

Parameters:
  • y (tuple) – Tuple with the suceptible and infected data.
  • t (array) – The time respective to each y set of samples.
  • Beta (float) – The Beta parameter.
  • r (float) – The r parameter.
  • sigma (float) – The sigma parameter.
Returns:

The derivative of the suceptible and infected data.

Return type:

tuple

models.differential_models.SIR(self, y, t, parameters, *args)

The function that computes the diferential set of equations of the SIR Epidemic Model.

Parameters:
  • y (tuple) – Tuple with the suceptible and infected data.
  • t (array) – The time respective to each y set of samples.
  • Beta (float) – The Beta parameter.
  • r (float) – The r parameter.
Returns:

The derivative of the suceptible and infected data.

Return type:

tuple

models.differential_models.SIRD(self, y, t, parameters)

The function that computes the diferential set of equations of the SIRD Epidemic Model.

Parameters:
  • y (tuple) – Tuple with the suceptible and infected data.
  • t (array) – The time respective to each y set of samples.
  • Beta (float) – The Beta parameter.
  • r (float) – The r parameter.
  • mi (float) – The mi parameter.
Returns:

The derivative of the suceptible and infected data.

Return type:

tuple

Module contents

copyright:2010 Marcelo Lima
license:BSD-3-Clause

Vanderlei Cunha Parro user image

LinkedIn Link Github Link

Gmail Link Youtube Link

Marcelo Mendes Lafetá Lima user image

LinkedIn Link Github Link Gmail Link

Felipe Brandão Ippolito user image

Github Link Gmail Link

Felipe Antonio Silva de Andrade user image

LinkedIn Link Github Link Gmail Link

Indices e tabelas