Моделирање Coronavirus податоци

Пред да започнеме со data fitting, треба да разгледаме кои параметри ги имаме и како новите равенки влијаат на кодот и како може да го поправиме.

Ова се сите параметри кои моментално му требаат на нашиот модел:

$N$: вкупна популација
$\beta$(t): очекуваниот број на луѓе на кои заразена личност ја пренесува болеста на ден
$D$: бројот на денови што личноста ја чува во себе и може да ја шири болеста
$\gamma$: дел од заразени личности што заздравуваат на ден ($\gamma= \frac{1}{D}$)
$R_{0_{start}}$, $x_0$, $R_{0_{end}}$, $k$ потребни за равенките на $R_0(t)$
$f$ - фактор параметар во $Beds(t)$
$Beds_0$ параметар во $R_0(t)$
$\delta$: времетраење на инкубациски период
$P(I\rightarrow C)$: веројатност за транзиција од состојба на инфициран во состојба на критичен
$P(C\rightarrow D)$: веројатност за починување додека личност е во критична здраствена состојба

Сега да разгледаме кои ни се потребни за fitting: Дефинитивно можеме да воочиме дека $N$: вкупна популација не ни треба да ја совпаѓаме бидејќи ја гледаме само популацијата на местото каде сакаме да го воспоставиме моделот и таа е константна, истото важи и за бројот на кревети за критични пациенти $Beds_0$. Факторите за транзиции $\delta, \gamma$ се фиксни вредности истотака: $\delta = 1/9$, $\gamma = 1/3$, кои се најверодостојните параметри кои ги најдовме низ примери од кодови за моделирање на слични сценарија [7]. Сларната вредност за факторот $f$ за бројот кревети не игра многу голема улога, повеќе го користиме бидејќи согледавме дека во примери служи за споредба како order of magnitude фактори влијаат на крајни сценарија, меѓутоа може и ќе учествува во fitting. Параметарот $\beta$, во последните равенки го бележиме со нотацијата $\beta (t)$ бидејќи очекуваниот број на луѓе на кои заразена личност ја пренесува болеста на ден ја пресметуваме преку $R_0 (t), \gamma$ бидејќи $R_0$ се менува не мора да имаме фиксен параметар за $\beta$ како по претходните тетратки.

Веројатностите $P(I\rightarrow C)$ и $P(C\rightarrow D)$, се поделени по возрасни групи (слично како $\alpha$ во оваа тетратка) што ни овозможува да правиме поделби врз основа на како наслението е распределено генереациски.

Што го намалува списокот на параметри за fitting на следните: $R_{0_start}$, $R_{0_end}$, $k$, $x_0$, $P(I\rightarrow C)$, $P(C\rightarrow D)$, $f$

Податоци за Coronavirus

Повеќето податоци што ќе ги користиме ги има низ повеќе примери на dashboards, податоци за возрасни групи, веројатности, кревети, $\alpha$ се земени од UNdata и од податоците за dashboard-от на Johs Hopkins - School of Engineering. [8]

Прост пример ова се најновите бројки за бројот на Починати, $D(t)$ во светот (овој график е правен со помош на Pandas функциите бидејќи нуди подобра опција за зумирање што сметаме дека е поубава за вакви примери):

import numpy as np
import pandas as pd
import lmfit
import warnings
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import mpld3
from lmfit.lineshapes import gaussian, lorentzian
from scipy.integrate import odeint
import plotly.graph_objects as go
import plotly.tools as tls
from plotly.offline import plot, iplot, init_notebook_mode
from IPython.core.display import display, HTML
init_notebook_mode(connected = True)
config={'showLink': False, 'displayModeBar': False}
mpld3.enable_notebook()
%matplotlib inline 
warnings.filterwarnings('ignore')

# read data
beds = pd.read_csv("https://raw.githubusercontent.com/hf2000510/infectious_disease_modelling/master/data/beds.csv", header=0)
agegroups = pd.read_csv("https://raw.githubusercontent.com/hf2000510/infectious_disease_modelling/master/data/agegroups.csv")
probabilities = pd.read_csv("https://raw.githubusercontent.com/hf2000510/infectious_disease_modelling/master/data/probabilities.csv")
covid_data = pd.read_csv("https://tinyurl.com/t59cgxn", parse_dates=["Date"], skiprows=[1])

beds_lookup = dict(zip(beds["Country"], beds["ICU_Beds"]))
agegroup_lookup = dict(zip(agegroups['Location'], agegroups[['0_9', '10_19', '20_29', '30_39', '40_49', '50_59', '60_69', '70_79', '80_89', '90_100']].values))

prob_I_to_C_1 = list(probabilities.prob_I_to_ICU_1.values)
prob_I_to_C_2 = list(probabilities.prob_I_to_ICU_2.values)
prob_C_to_Death_1 = list(probabilities.prob_ICU_to_Death_1.values)
prob_C_to_Death_2 = list(probabilities.prob_ICU_to_Death_2.values)

plt.clf
covid_data.rename(columns = {'Value':'Починати'}, inplace = True)
covid_data.groupby("Date").sum()[["Починати"]].plot(
    figsize=(6, 3), 
    title="Covid-19 вкупен број починати (светот)",
    color="#696969"
)
covid_data.rename(columns = {'Починати':'Value'}, inplace = True)

Графикот кои го користиме е пример што го најдовме на интернет, каде користиме всушност само податоци (dataset) за бројот на починати а не на бројот на пријавени. Причината е веројатно бидејќи шум кој би се јавил во податоците поврзан со пријавените случаеви и бројот на тестови (без разлика што има голем број на тестови еве и најголем да е, самите држави не ги тестираат сите кои се инфицирани). Зошто сметаме дека без разлика колку е добар моделот, надворешните фактори играат многу голема улога [9]. Пример сме држава која правела по 200 теста на ден, доколку нагло добиеме доволно тестови да ги тестираме сите луѓе не значи дека следниот ден ќе имаме 50% повеќе заболени, туку имплицира дека бројот на тестови е зголемен, што од оваа перспектива на луѓе што се запознааваат со овој свет изгледа фантастично како воопшто луѓе даваат добри предвидувања. Да не речам и невозможно. Секој случај, бројот на починати истотака може да биде под влијание од системите на државите но бројките што кружат низ интернетот се најверојатно блиску до вистинските.

Програмирање на проширениот модел

Имајќи ги во предвид равенките погоре наведни, сега може да го програмираме моделот слично како претходните тетратки. Ова е резултатот од сценариото каде нема доволно кревети за критични случаеви:

import numpy as np
import pandas as pd
import lmfit
import warnings
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import mpld3
from lmfit.lineshapes import gaussian, lorentzian
from scipy.integrate import odeint
import plotly.graph_objects as go
import plotly.tools as tls
from plotly.offline import plot, iplot, init_notebook_mode
from IPython.core.display import display, HTML
init_notebook_mode(connected = True)
config={'showLink': False, 'displayModeBar': False}
mpld3.enable_notebook()
%matplotlib inline 
warnings.filterwarnings('ignore')

def deriv(y, t, beta, gamma, sigma, N, p_I_to_C, p_C_to_D, Beds):
    S, E, I, C, R, D = y

    dSdt = -beta(t) * I * S / N
    dEdt = beta(t) * I * S / N - sigma * E
    dIdt = sigma * E - 1/12.0 * p_I_to_C * I - gamma * (1 - p_I_to_C) * I
    dCdt = 1/12.0 * p_I_to_C * I - 1/7.5 * p_C_to_D * min(Beds(t), C) - max(0, C-Beds(t)) - (1 - p_C_to_D) * 1/6.5 * min(Beds(t), C)
    dRdt = gamma * (1 - p_I_to_C) * I + (1 - p_C_to_D) * 1/6.5 * min(Beds(t), C)
    dDdt = 1/7.5 * p_C_to_D * min(Beds(t), C) + max(0, C-Beds(t))
    return dSdt, dEdt, dIdt, dCdt, dRdt, dDdt

def logistic_R_0(t, R_0_start, k, x0, R_0_end):
    return (R_0_start-R_0_end) / (1 + np.exp(-k*(-t+x0))) + R_0_end

def beta(t):
        return logistic_R_0(t, R_0_start, k, x0, R_0_end) * gamma

    
def Beds(t):
    beds_0 = beds_per_100k / 100_000 * N
    return beds_0 + s*beds_0*t  # 0.003

def logistic_R_0(t, R_0_start, k, x0, R_0_end):
    return (R_0_start-R_0_end) / (1 + np.exp(-k*(-t+x0))) + R_0_end

gamma = 1.0/9.0
sigma = 1.0/3.0

days=500
agegroups=[100000, 100000, 100000, 100000, 100000, 100000, 100000, 100000, 100000]
beds_per_100k=50
R_0_start=4.0
k=1.0
x0=60
R_0_end=1.0 
prob_I_to_C=0.05
prob_C_to_D=0.6
s=0.003

N = sum(agegroups)

y0 = N-1.0, 1.0, 0.0, 0.0, 0.0, 0.0
t = np.linspace(0, days-1, days)
ret = odeint(deriv, y0, t, args=(beta, gamma, sigma, N, prob_I_to_C, prob_C_to_D, Beds))
S, E, I, C, R, D = ret.T
R_0_over_time = [beta(i)/gamma for i in range(len(t))]

data_com = []

data_1 = go.Scatter(x = t, 
                  y = S, 
                  mode = 'lines',
                  visible = True,
                  line=dict(color="blue",
                            width=2),
                  name = "Подлежни; <i>S(t)</i>",
                  hovertemplate = '<i> %{y:.0f} </i> луѓе')


data_2 = go.Scatter(x = t, 
                  y = E, 
                  mode = 'lines',
                  visible = True,
                  line=dict(color="#CCCC00",
                            width=2),
                  name = "Изложени; <i>E(t)</i>",
                  hovertemplate = '<i> %{y:.0f} </i> луѓе')


data_3 = go.Scatter(x = t, 
                  y = I,
                  mode = 'lines',
                  visible = True, 
                  line=dict(color="red",
                            width=2),
                  name = "Заразени; <i>I(t)</i>",
                  hovertemplate = '<i> %{y:.0f} </i> луѓе')


data_4 = go.Scatter(x = t,
                  y = R, 
                  mode = 'lines',
                  visible = True,
                  line=dict(color="green",
                            width=2),
                  name = "Оздравени; <i>R(t)</i>",
                  hovertemplate = '<i> %{y:.0f} </i> луѓе')


data_5 = go.Scatter(x = t,
                  y =  D, 
                  visible = True, 
                  mode = 'lines',
                  line=dict(color="#696969",
                            width=2),
                  name = "Починати; <i>D(t)</i>",
                  hovertemplate = '<i> %{y:.0f} </i> луѓе')

data_6 = [go.Scatter(x = t,
                  y =  C, 
                  visible = True, 
                  mode = 'lines',
                  line=dict(color="tomato",
                            width=2),
                  name = "Критични; <i>C(t)</i>",
                  hovertemplate = '<i> %{y:.0f} </i> луѓе')]

data_com.append(data_1)
data_com.append(data_2)
data_com.append(data_3)
data_com.append(data_4)
data_com.append(data_5)

data_R0 = [dict(visible = True,
                    x = t,
                    y = R_0_over_time,
                    name = 'R<sub>0</sub>',
                    hoverinfo = "all",
                    hovertemplate = "%{y: .1f}",
                    line=dict(width=1.5,
                              dash = 'dot',
                              color = 'rgb(0, 204, 150)'),
                    xaxis='x2',
                    yaxis='y2',
                    showlegend = False
               )]

total_CFR = [0] + [100 * D[i] / sum(sigma*E[:i]) if sum(sigma*E[:i])>0 else 0 for i in range(1, len(t))]
daily_CFR = [0] + [100 * ((D[i]-D[i-1]) / ((R[i]-R[i-1]) + (D[i]-D[i-1]))) if max((R[i]-R[i-1]), (D[i]-D[i-1]))>10 else 0 for i in range(1, len(t))]

data_alpha1 = [dict(visible = False,
                    x = t,
                    y = total_CFR,
                    name = 'Вкупно &#120572;',
                    hoverinfo = "all",
                    hovertemplate = "%{y: .2f}",
                    line=dict(width=1.5,
                              dash = 'dot',
                              color = 'rgb(238,144,238)'),
                    xaxis='x2',
                    yaxis='y2',
                    showlegend = False,
               )]

data_alpha2 = [dict(visible = False,
                    x = t,
                    y = daily_CFR,
                    name = 'Дневно &#120572; ',
                    hoverinfo = "all",
                    hovertemplate = "%{y: .2f}",
                    line=dict(width=1.5,
                              dash = 'dot',
                              color = 'palegreen'),
                    xaxis='x2',
                    yaxis='y2',
                    showlegend = False,
               )]

newDs = [0] + [D[i]-D[i-1] for i in range(1, len(t))]

data_death1 = [dict(visible = False,
                    x = t,
                    y = newDs,
                    name = 'Вкупно починати',
                    hoverinfo = "all",
                    hovertemplate = "%{y}",
                    line=dict(width=1.5,
                              dash = 'dot',
                              color = 'gray'),
                    xaxis='x2',
                    yaxis='y2',
                    showlegend = False,
               )]

data_death2 = [dict(visible = False,
                    x = t,
                    y = [max(0, C[i]-Beds(i)) for i in range(len(t))],
                    name = 'Починати од недостиг ресурси',
                    hoverinfo = "all",
                    hovertemplate = "%{y}",
                    line=dict(width=1.5,
                              dash = 'dot',
                              color = 'tomato'),
                    xaxis='x2',
                    yaxis='y2',
                    showlegend = False,
               )]

    
data = data_com + data_R0 + data_alpha1 + data_alpha2  + data_death1 + data_death2 + data_6

# Setup the layout of the figure
layout = go.Layout(
    updatemenus=[
        dict(
            active = 0, 
            x=0.2,
            y=1.1,
            yanchor="top",
            buttons=list([
                dict(label="(1) R<sub>0</sub> преку време",
                             method="update",
                             args=[{"visible": [True, True, True, True, True, True, False, False, False, False, True]},
                                   {'annotations': [
                                               dict(text="Прикажи: ", 
                                                     showarrow=False,
                                                     x=0.25, 
                                                     y=1.15, 
                                                     yref="paper", 
                                                     align="left"),
                                               dict(x=0.15,
                                                    y=0.76,
                                                    showarrow=False,
                                                    text='R<sub>0</sub> преку време',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper'),
                                                dict(x=0.228,
                                                    y=0.21,
                                                    showarrow=False,
                                                    text='t',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper'),
                                                dict(x=0.01,
                                                    y=0.45,
                                                    showarrow=False,
                                                    textangle = -90,
                                                    text='R<sub>0</sub>',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper')] ,
                                    'xaxis2': dict(range = [0,500],
                                                   tickvals = [0, 100, 200, 300, 400, 500],
                                                   linecolor='#000',
                                                   domain=[0.08, 0.38],
                                                   anchor='y2',
                                                   mirror = True,
                                                   side='bottom',
                                                   ticks='',
                                                   showline=True),
                                   'yaxis2': dict(autorange=True,
                                                  range=[0.5, 4.5],
                                                  tickvals = [1, 2, 3, 4],
                                                  linecolor='#000',
                                                  domain=[0.30, 0.70],
                                                  anchor='x2',
                                                  mirror = True,
                                                  ticks='',
                                                  showline=True)}
                                  ]) ,
                
                dict(label="(2) &#120572; преку време",
                             method="update",
                             args=[{"visible": [True, True, True, True, True, False, True, True, False, False, True]},
                                   {'annotations': [
                                               dict(text="Прикажи: ", 
                                                     showarrow=False,
                                                     x=0.25, 
                                                     y=1.15, 
                                                     yref="paper", 
                                                     align="left"),
                                               dict(x=0.15,
                                                    y=0.76,
                                                    showarrow=False,
                                                    text='&#120572; преку време',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper'),
                                               dict(x=0.228,
                                                    y=0.21,
                                                    showarrow=False,
                                                    text='t',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper'),
                                               dict(x=0.013,
                                                    y=0.48,
                                                    showarrow=False,
                                                    textangle = -90,
                                                    text='&#120572;',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper')],
                                    'xaxis2': dict(range = [0,500],
                                                   tickvals = [0, 100, 200, 300, 400, 500],
                                                   linecolor='#000',
                                                   domain=[0.08, 0.38],
                                                   anchor='y2',
                                                   mirror = True,
                                                   side='bottom',
                                                   ticks='',
                                                   showline=True,
                                                   tickfont = dict(size=9)),
                                   'yaxis2': dict(autorange = False,
                                                  range = [0, 2.8],
                                                  tickvals = [0, 0.4,  0.8,  1.2,  1.6,  2.0,  2.4, 2.8],
                                                  linecolor='#000',
                                                  domain=[0.30, 0.70],
                                                  anchor='x2',
                                                  mirror = True,
                                                  ticks='',
                                                  showline=True)}
                                  ]) , 
                dict(label="(3) Смртни случаеви",
                             method="update",
                             args=[{"visible": [True, True, True, True, True, False, False, False, True, True, True]},
                                   {'annotations': [
                                               dict(text="Прикажи: ", 
                                                     showarrow=False,
                                                     x=0.25, 
                                                     y=1.15, 
                                                     yref="paper", 
                                                     align="left"),
                                               dict(x=0.125,
                                                    y=0.76,
                                                    showarrow=False,
                                                    text='Починати дневно',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper'),
                                               dict(x=0.228,
                                                    y=0.21,
                                                    showarrow=False,
                                                    text='t',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper'),
                                               dict(x=0.017,
                                                    y=0.48,
                                                    showarrow=False,
                                                    textangle = -90,
                                                    text='#',
                                                    font=dict(size=10),
                                                    xref='paper',
                                                    yref='paper')],
                                    'xaxis2': dict(range = [0,500],
                                                   tickvals = [0, 100, 200, 300, 400, 500],
                                                   linecolor='#000',
                                                   domain=[0.08, 0.38],
                                                   anchor='y2',
                                                   mirror = True,
                                                   side='bottom',
                                                   ticks='',
                                                   showline=True,
                                                  tickfont = dict(size=9)),
                                   'yaxis2': dict(autorange = False,
                                                  range = [0, 75],
                                                  tickvals = [0, 10, 20, 30, 40, 50, 60, 70],
                                                  linecolor='#000',
                                                  domain=[0.30, 0.70],
                                                  anchor='x2',
                                                  mirror = True,
                                                  ticks='',
                                                  showline=True)}
                                  ]) 
                ]),
            direction="down"
            )],
                  title = "Проширен SIR модел",
                  title_x = 0.5, 
                  xaxis_title='t',
                  annotations=[
                               dict(text="Прикажи: ", 
                                     showarrow=False,
                                     x=0.25, 
                                     y=1.15, 
                                     yref="paper", 
                                     align="left"),
                               dict(x=0.15,
                                    y=0.76,
                                    showarrow=False,
                                    text='R<sub>0</sub> преку време',
                                    font=dict(size=10),
                                    xref='paper',
                                    yref='paper'),
                                dict(x=0.228,
                                    y=0.21,
                                    showarrow=False,
                                    text='t',
                                    font=dict(size=10),
                                    xref='paper',
                                    yref='paper'),
                                dict(x=0.027,
                                    y=0.48,
                                    showarrow=False,
                                    textangle = -90,
                                    text='R<sub>0</sub>',
                                    font=dict(size=10),
                                    xref='paper',
                                    yref='paper')],
    
                  xaxis=dict(range=[0,500], 
                             mirror=False,
                             ticks='outside',
                             showline=True,
                             tickvals = [0, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500],
                             linecolor='#000',
                             tickfont = dict(size=11)),
                  yaxis_title='Број на луѓе',
                  yaxis=dict(range=[0,940000], 
                             mirror=False,
                             ticks='outside', 
                             showline=True,
                             showspikes = True,
                             linecolor='#000',
                             tickvals = [0, 100000, 200000, 300000, 400000, 500000, 600000, 700000, 800000, 900000], 
                             tickfont = dict(size=11)),
                  xaxis2=dict(range = [0,500],
                              tickvals = [0, 100, 200, 300, 400, 500],
                              linecolor='#000',
                                domain=[0.09, 0.38],
                                anchor='y2',
                                mirror = True,
                                side='bottom',
                                ticks='',
                                showline=True),
                  yaxis2=dict(autorange=False,
                              range = [0.5, 4.5],
                              tickvals = [1, 2, 3, 4],
                            linecolor='#000',
                            domain=[0.30, 0.70],
                            anchor='x2',
                            mirror = True,
                            ticks='',
                            showline=True),
                  plot_bgcolor='#fff', 
                  hovermode = 'x unified',
                  width = 700, 
                  height = 400,
                  font = dict(size = 10),
                  margin=go.layout.Margin(l=50,
                                         r=50,
                                         b=60,
                                         t=35))

# Plot function saves as html or with ipplot
fig_extend = dict(data=data, layout=layout)
plot(fig_extend, filename = 'fig_extend.html', config = config)
display(HTML('fig_extend.html'))

На третиот мал график може да се истражува како бројот на смртни случаеви се менува дневно како резултат на недостигот на кревети. Сега останува уште да направиме Curve Fitting за Италија во следната тетратка.