Usage: SIR-derived models

Open In Colab

Here, we will create example datasets with simulated values of SIR-derived models. Then, we will perform scenario analysis with them.

Preparation

Prepare tha packages and create a instance of PopulationData class to save population values.

[1]:
# !pip install covsirphy --upgrade
from pprint import pprint
import covsirphy as cs
cs.__version__
[1]:
'2.20.3-gamma'

Create example dataset with theoretical values

We will use ExampleData class to perform simulation with preset initial values and parameters. \(\tau\) (coeficient for non-dimensionalization) will be set as \(1440\ \mathrm{[min]}\). The first date of records will be 01Jan2020 as an example.

[2]:
# Set tau value and start date of records
example_data = cs.ExampleData(tau=1440, start_date="01Jan2020")

No records were registered at this time.

[3]:
# Check records
example_data.cleaned()
[3]:
Date Country Province Confirmed Infected Fatal Recovered

ExampleData class is a child class of JHUData. i.e. We can use the example data in scenario analysis. Example codes will be shown in “Scenario analysis with theoretical data” subsection.

[4]:
issubclass(cs.ExampleData, cs.JHUData)
[4]:
True
[5]:
isinstance(example_data, cs.JHUData)
[5]:
True

SIR model

Let’s start with the simlest SIR model. “Susceptible people” may meet “Infected” persons and may be confirmed as “Infected”. “Infected” patients will move to “Recovered” compertment later.

\begin{align*} \mathrm{S} \overset{\beta I}{\longrightarrow} \mathrm{I} \overset{\gamma}{\longrightarrow} \mathrm{R} \\ \end{align*}

Variables:

  • \(\mathrm{S}\): Susceptible (= Population - Confirmed)

  • \(\mathrm{I}\): Infected (=Confirmed - Recovered - Fatal)

  • \(\mathrm{R}\): Recovered or Fatal (= Recovered + Fatal)

Parameters:

  • \(\beta\): Effective contact rate \(\mathrm{[1/min]}\)

  • \(\gamma\): Recovery (+ Mortality) rate \(\mathrm{[1/min]}\)

Note:
Though \(R\) in SIR model is “Recovered and have immunity”, we defines \(R\) as “Recovered or fatal”. This is because mortality rate cannot be ignored in our COVID-19 outbreak.

Non-dimensional SIR model

To simplify the model, we will remove the units of the variables from the ODE model.

Set \((S, I, R) = N \times (x, y, z)\) and \((T, \beta, \gamma) = (\tau t, \tau^{-1} \rho, \tau^{-1} \sigma)\).

This results in the ODE

\begin{align*} & \frac{\mathrm{d}x}{\mathrm{d}t}= - \rho x y \\ & \frac{\mathrm{d}y}{\mathrm{d}t}= \rho x y - \sigma y \\ & \frac{\mathrm{d}z}{\mathrm{d}t}= \sigma y \\ \end{align*}

Where \(N\) is the total population and \(\tau\) is a coefficient ([min], is an integer to simplify).

The model name and preset of parameter and initial values are registed as class variables of SIR class.

[6]:
# Model name
print(cs.SIR.NAME)
# Example parameter values
pprint(cs.SIR.EXAMPLE, compact=True)
SIR
{'param_dict': {'rho': 0.2, 'sigma': 0.075},
 'population': 1000000,
 'step_n': 180,
 'y0_dict': {'Fatal or Recovered': 0, 'Infected': 1000, 'Susceptible': 999000}}

With the preset values, ExampleData instance will produce a example data.

[7]:
model = cs.SIR
area = {"country": "Full", "province": model.NAME}
# Add records with SIR model
example_data.add(model, **area)

We can get example records with ExampleData.specialized() method.

[8]:
# Records with model variables
df = example_data.specialized(model, **area)
df.head()
[8]:
Date Susceptible Infected Fatal or Recovered
0 2020-01-02 998787 1133 80
1 2020-01-03 998547 1283 170
2 2020-01-04 998273 1454 273
3 2020-01-05 997964 1647 389
4 2020-01-06 997614 1865 521

With covsirphy.line_plot() function, figures will be shown (or saved when filename argument was applied).

[9]:
# Line plot with the example data
cs.line_plot(df.set_index("Date"), title=f"Example data of {model.NAME} model", y_integer=True)
_images/usage_theoretical_19_0.png

Reproduction number

Reproduction number of SIR model is defined as follows.

\begin{align*} R_0 = \rho \sigma^{-1} = \beta \gamma^{-1} \end{align*}

\(R_0\) (“R naught”) means “the average number of secondary infections caused by an infected host” (external link: Infection Modeling — Part 1). When \(x=\frac{1}{R_0}\), \(\frac{\mathrm{d}y}{\mathrm{d}t}=0\) (the number of infected cases does not change).

We can calculate reproduction number using .calc_r0() method.

[10]:
# Calculate reproduction number
# Note: population value will be applied, but not used in calculation
param_dict = cs.SIR.EXAMPLE["param_dict"].copy()
model_instance = cs.SIR(population=100000, **param_dict)
r0 = model_instance.calc_r0()
print(f"Reproduction number of {model_instance.NAME} model: {r0}")
Reproduction number of SIR model: 2.67

SIR-D model

Because we are measuring the number of fatal cases and recovered cases separately, we can use two variables (“Recovered” and “Deaths”) instead of “Recovered + Deaths” in the mathematical model. We call this model as SIR-D model.

\begin{align*} \mathrm{S} \overset{\beta I}{\longrightarrow}\ & \mathrm{I} \overset{\gamma}{\longrightarrow} \mathrm{R} \\ & \mathrm{I} \overset{\alpha}{\longrightarrow} \mathrm{D} \\ \end{align*}

Variables:

  • \(\mathrm{S}\): Susceptible (= Population - Confirmed)

  • \(\mathrm{I}\): Infected (=Confirmed - Recovered - Fatal)

  • \(\mathrm{R}\): Recovered

  • \(\mathrm{D}\): Fatal

Parameters:

  • \(\alpha\): Mortality rate \(\mathrm{[1/min]}\)

  • \(\beta\): Effective contact rate \(\mathrm{[1/min]}\)

  • \(\gamma\): Recovery rate \(\mathrm{[1/min]}\)

Non-dimensional SIR-D model

Set \((S, I, R, D) = N \times (x, y, z, w)\) and \((T, \alpha, \beta, \gamma) = (\tau t, \tau^{-1} \kappa, \tau^{-1} \rho, \tau^{-1} \sigma)\). This results in the ODE

\begin{align*} & \frac{\mathrm{d}x}{\mathrm{d}t}= - \rho x y \\ & \frac{\mathrm{d}y}{\mathrm{d}t}= \rho x y - (\sigma + \kappa) y \\ & \frac{\mathrm{d}z}{\mathrm{d}t}= \sigma y \\ & \frac{\mathrm{d}w}{\mathrm{d}t}= \kappa y \\ \end{align*}

The model name and preset values are registered in SIRD class.

[11]:
# Model name
print(cs.SIRD.NAME)
# Example parameter values
pprint(cs.SIRD.EXAMPLE, compact=True)
SIR-D
{'param_dict': {'kappa': 0.005, 'rho': 0.2, 'sigma': 0.075},
 'population': 1000000,
 'step_n': 180,
 'y0_dict': {'Fatal': 0,
             'Infected': 1000,
             'Recovered': 0,
             'Susceptible': 999000}}

Example data is here.

[12]:
model = cs.SIRD
area = {"country": "Full", "province": model.NAME}
# Add records with SIR model
example_data.add(model, **area)
# Records with model variables
df = example_data.specialized(model, **area)
cs.line_plot(df.set_index("Date"), title=f"Example data of {model.NAME} model", y_integer=True)
_images/usage_theoretical_26_0.png

Reproduction number

Reproduction number of SIR-D model is defined as follows.

\begin{align*} R_0 = \rho (\sigma + \kappa)^{-1} = \beta (\gamma + \alpha)^{-1} \end{align*}

We can calculate reproduction number using .calc_r0() method.

[13]:
# Calculate reproduction number
# Note: population value will be applied, but not used in calculation
param_dict = cs.SIRD.EXAMPLE["param_dict"].copy()
model_instance = cs.SIRD(population=100000, **param_dict)
r0 = model_instance.calc_r0()
print(f"Reproduction number of {model_instance.NAME} model: {r0}")
Reproduction number of SIR-D model: 2.5

SIR-F model

In the initial phase of COVID-19 outbreak, many cases were confirmed after they died. To consider this issue, “S + I \(\to\) Fatal + I” should be added. We call the next model as SIR-F model. When \(\alpha_{1}=0\), no difference with the SIR-D model.

\begin{align*} \mathrm{S} \overset{\beta I}{\longrightarrow} \mathrm{S}^\ast \overset{\alpha_1}{\longrightarrow}\ & \mathrm{F} \\ \mathrm{S}^\ast \overset{1 - \alpha_1}{\longrightarrow}\ & \mathrm{I} \overset{\gamma}{\longrightarrow} \mathrm{R} \\ & \mathrm{I} \overset{\alpha_2}{\longrightarrow} \mathrm{F} \\ \end{align*}

Variables:

  • \(\mathrm{S}\): Susceptible (= Population - Confirmed)

  • \(\mathrm{S}^\ast\): Confirmed and un-categorized

  • \(\mathrm{I}\): Confirmed and categorized as Infected

  • \(\mathrm{R}\): Confirmed and categorized as Recovered

  • \(\mathrm{F}\): Confirmed and categorzied as Fatal

Parameters:

  • \(\alpha_1\): Direct fatality probability of \(\mathrm{S}^\ast\) (non-dimensional)

  • \(\alpha_2\): Mortality rate of Infected cases \(\mathrm{[1/min]}\)

  • \(\beta\): Effective contact rate \(\mathrm{[1/min]}\)

  • \(\gamma\): Recovery rate \(\mathrm{[1/min]}\)

Notes on \(\mathrm{S}^\ast\) variable:
\(\mathrm{S}^\ast\) describes the cases who are actually carriers of the disease without anyone (including themselves) knowing about it, who either die and they are confirmed positive after death, while some others are moved to infected after being confirmed.

In JHU-style dataset, we know the number of cases who were confirmed with COVID-19, but we do not know the number of died cases who died without COVID-19. Essentially \(\mathrm{S}^\ast\) serves as an auxiliary compartment in SIR-F model to separate the two death situations and insert a probability factor of {\(\alpha_1\), \(1 - \alpha_1\)}.

Notes on the difference of SIR-D and SIR-F model:
\(\alpha_1\) is small at this time because performance of PCR tests was improved, but we can use SIR-F model rather than SIR-D model as an enhanced model even now becase \(\alpha_1\) can be 0 in the ODE model.

SIR-F model was developed with Kaggle: COVID-19 data with SIR model.

Non-dimensional SIR-F model

Set \((S, I, R, F) = N \times (x, y, z, w)\) and \((T, \alpha_1, \alpha_2, \beta, \gamma) = (\tau t, \theta, \tau^{-1} \kappa, \tau^{-1} \rho, \tau^{-1} \sigma)\). This results in the ODE

\begin{align*} & \frac{\mathrm{d}x}{\mathrm{d}t}= - \rho x y \\ & \frac{\mathrm{d}y}{\mathrm{d}t}= \rho (1-\theta) x y - (\sigma + \kappa) y \\ & \frac{\mathrm{d}z}{\mathrm{d}t}= \sigma y \\ & \frac{\mathrm{d}w}{\mathrm{d}t}= \rho \theta x y + \kappa y \\ \end{align*}

We use SIRF class.

[14]:
# Model name
print(cs.SIRF.NAME)
# Example parameter values
pprint(cs.SIRF.EXAMPLE, compact=True)
SIR-F
{'param_dict': {'kappa': 0.005, 'rho': 0.2, 'sigma': 0.075, 'theta': 0.002},
 'population': 1000000,
 'step_n': 180,
 'y0_dict': {'Fatal': 0,
             'Infected': 1000,
             'Recovered': 0,
             'Susceptible': 999000}}

Example data is here.

[15]:
model = cs.SIRF
area = {"country": "Full", "province": model.NAME}
# Add records with SIR model
example_data.add(model, **area)
# Records with model variables
df = example_data.specialized(model, **area)
cs.line_plot(df.set_index("Date"), title=f"Example data of {model.NAME} model", y_integer=True)
_images/usage_theoretical_36_0.png

Reproduction number

Reproduction number of SIR-F model is defined as follows.

\begin{align*} R_0 = \rho (1 - \theta) (\sigma + \kappa)^{-1} = \beta (1 - \alpha_1) (\gamma + \alpha_2)^{-1} \end{align*}

We can calculate reproduction number using .calc_r0() method.

[16]:
# Calculate reproduction number
# Note: population value will be applied, but not used in calculation
param_dict = cs.SIRF.EXAMPLE["param_dict"].copy()
model_instance = cs.SIRF(population=100000, **param_dict)
r0 = model_instance.calc_r0()
print(f"Reproduction number of {model_instance.NAME} model: {r0}")
Reproduction number of SIR-F model: 2.5

SIR-F with exposed/waiting cases

The next model is SEWIR-F model. The number of exposed cases in latent period (E) and wating cases for confirmation (W) are un-measurable variables, but key variables as well as S, I, R, F. If E and W are large, outbreak will occur in the near future. Let’s replace S\(\overset{\beta I}{\longrightarrow}\)S\(^\ast\) as follows because W also has infectivity.

\begin{align*} \mathrm{S} \overset{\beta_1 (W+I)}{\longrightarrow} \mathrm{E} \overset{\beta_2}{\longrightarrow} \mathrm{W} \overset{\beta_3}{\longrightarrow} \mathrm{S}^\ast \overset{\alpha_1}{\longrightarrow}\ & \mathrm{F} \\ \mathrm{S}^\ast \overset{1 - \alpha_1}{\longrightarrow}\ & \mathrm{I} \overset{\gamma}{\longrightarrow} \mathrm{R} \\ & \mathrm{I} \overset{\alpha_2}{\longrightarrow} \mathrm{F} \\ \end{align*}

Variables:

  • \(\mathrm{S}\): Susceptible

  • \(\mathrm{E}\): Exposed and in latent period (without infectivity)

  • \(\mathrm{W}\): Waiting for confirmaion fiagnosis (with infectivity)

  • \(\mathrm{S}^\ast\): Confirmed and un-categorized

  • \(\mathrm{I}\): Confirmed and categorized as Infected

  • \(\mathrm{R}\): Confirmed and categorized as Recovered

  • \(\mathrm{F}\): Confirmed and categorzied as Fatal

Parameters:

  • \(\alpha_1\): Direct fatality probability of \(\mathrm{S}^\ast\) (non-dimensional)

  • \(\alpha_2\): Mortality rate of Infected cases \(\mathrm{[1/min]}\)

  • \(\beta_1\): Exposure rate (the nymber of encounter with the virus in a minute) \(\mathrm{[1/min]}\)

  • \(\beta_2\): Inverse of latent period \(\mathrm{[1/min]}\)

  • \(\beta_3\): Inverse of waiting time for confirmation \(\mathrm{[1/min]}\)

  • \(\gamma\): Recovery rate \(\mathrm{[1/min]}\)

Non-dimensional SEWIR-F model

Set \((S, E, W, I, R, F) = N \times (x_1, x_2, x_3, y, z, w)\), \((T, \alpha_1) = (\tau t, \theta)\) and \((\alpha_2, \beta_i, \gamma) = \tau^{-1} \times (\kappa, \rho_i, \sigma)\). This results in the ODE

\begin{align*} & \frac{\mathrm{d}x_1}{\mathrm{d}t}= - \rho_1 x_1 (x_3 + y) \\ & \frac{\mathrm{d}x_2}{\mathrm{d}t}= \rho_1 x_1 (x_3 + y) - \rho_2 x_2 \\ & \frac{\mathrm{d}x_3}{\mathrm{d}t}= \rho_2 x_2 - \rho_3 x_3 \\ & \frac{\mathrm{d}y}{\mathrm{d}t}= (1-\theta) \rho_3 x_3 - (\sigma + \kappa) y \\ & \frac{\mathrm{d}z}{\mathrm{d}t}= \sigma y \\ & \frac{\mathrm{d}w}{\mathrm{d}t}= \theta \rho_3 x_3 + \kappa y \\ \end{align*}
Note:
We cannot use SEWIR-F model for parameter estimation because we do not have records of Exposed and Waiting. Please use SIR-F model with covsirphy.SIRF class.

SEWIRF class is for the SEWIR-F model.

[17]:
# Model name
print(cs.SEWIRF.NAME)
# Example parameter values
pprint(cs.SEWIRF.EXAMPLE, compact=True)
SEWIR-F
{'param_dict': {'kappa': 0.005,
                'rho1': 0.2,
                'rho2': 0.167,
                'rho3': 0.167,
                'sigma': 0.075,
                'theta': 0.002},
 'population': 1000000,
 'step_n': 180,
 'y0_dict': {'Exposed': 3000,
             'Fatal': 0,
             'Infected': 1000,
             'Recovered': 0,
             'Susceptible': 994000,
             'Waiting': 2000}}

Example records are here.

[18]:
model = cs.SEWIRF
area = {"country": "Full", "province": model.NAME}
# Add records with SIR model
example_data.add(model, **area)
# Records with model variables
df = example_data.specialized(model, **area)
cs.line_plot(df.set_index("Date"), title=f"Example data of {model.NAME} model", y_integer=True)
_images/usage_theoretical_45_0.png

Reproduction number

Reproduction number of SEWIR-F model is defined as follows.

\begin{align*} R_0 = \rho_1 /\rho_2 * \rho_3 (1-\theta) (\sigma + \kappa)^{-1} \end{align*}

We can calculate reproduction number using .calc_r0() method.

[19]:
# Calculate reproduction number
# Note: population value will be applied, but not used in calculation
param_dict = cs.SEWIRF.EXAMPLE["param_dict"].copy()
model_instance = cs.SEWIRF(population=100000, **param_dict)
r0 = model_instance.calc_r0()
print(f"Reproduction number of {model_instance.NAME} model: {r0}")
Reproduction number of SEWIR-F model: 2.49

SIR-F with vaccination

Vaccination is a key factor to prevent outbreak as you know.

In the previous version, we defined SIR-FV model with \(\omega\) (vaccination rate) and

\[\begin{split}\frac{\mathrm{d}S}{\mathrm{d}T}= - \beta S I - \omega N \\\end{split}\]

However, SIR-FV model was removed because vaccinated persons may move to the other compartments, including “Susceprtible”. Please use SIR-F model for simulation and parameter estimation with adjusted parameter values, considering the impact of vaccinations on infectivity, its effectivity and safety.

SIR-F with re-infection

Re-infection (Recovered -> Susceptble) is sometimes reported and we can consider SIR-S (SIR-FS) model. However, this is not impremented at this time because we do not have data regarding re-infection. SIR-F model could be the final model in our data-driven approach at this time.

Re-infection changes the parameter values of SIR-F model. There are two patterns.

  1. If re-infected case are counted as new confirmed cases and removed from “Recovered” compartment, \(\sigma\) will be decreased.

  2. If re-infected cases are counted as new confirmed cases and NOT removed from “Recovered” compartment, \(\rho\) will be increased because “Susceptible” will be decreased.

Impact of parameter change

Because ExampleData class is a subclass of JHUData, we can perform scenario analysis with example datasets easily. We evaluate the impact of parameter changes.

Here, we will use the following scenarios. For explanation, \(\tau=1440\), the start date is 01Jan2020, population is 1,000,000 and country name is “Theoretical”. Their scenarios are not based on actual data.

name

01Jan2020 - 31Jan2020

01Feb2020 - 31Dec2020

Main

SIR-F

SIR-F

Lockdown

SIR-F

SIR-F with 50% of \(\rho\)

Medicine

SIR-F

SIR-F with 50% of \(\kappa\) and 200% of \(\sigma\)

Vaccine

SIR-F

SIR-F with 80% of \(\rho\), 60% of \(\kappa\) and 120% of \(\sigma\)

As baseline (main scenario), we use preset values of the SIR-F model.

[20]:
# Preset of SIR-F parameters and initial values
preset_dict = cs.SIRF.EXAMPLE["param_dict"]
preset_dict
[20]:
{'theta': 0.002, 'kappa': 0.005, 'rho': 0.2, 'sigma': 0.075}

Create records from 01Jan2020 to 31Jan2020. These records will be used commonly in the scenarios.

[21]:
area = {"country": "Theoretical"}
# Create dataset from 01Jan2020 to 31Jan2020
example_data.add(cs.SIRF, step_n=30, **area)

Create Scenario instance for scenario analysis.

[22]:
# Create Scenario instance
snl = cs.Scenario(tau=1440, **area)
snl.register(example_data)

Then, confirm the records with Scenario.records() instance.

[23]:
# Show records with Scenario instance
record_df = snl.records()
display(record_df.head())
display(record_df.tail())
_images/usage_theoretical_58_0.png
Date Infected Fatal Recovered
0 2020-01-02 1127 6 80
1 2020-01-03 1270 12 169
2 2020-01-04 1430 19 271
3 2020-01-05 1612 28 385
4 2020-01-06 1816 37 513
Date Infected Fatal Recovered
25 2020-01-27 21210 940 13068
26 2020-01-28 23730 1061 14753
27 2020-01-29 26524 1196 16637
28 2020-01-30 29616 1347 18741
29 2020-01-31 33030 1515 21088
Note:
Record on 01Jan2020 was removed because the number of recovered cases is 0 and this sometimes causes error in estimation.

Then, set the records from 02Jan2020 to 31Jan2020 as the 0th phase. The 0th phase is commonly used in the all scenarios.

[24]:
# Set 0th phase from 02Jan2020 to 31Jan2020 with preset parameter values
snl.clear(include_past=True)
snl.add(end_date="31Jan2020", model=cs.SIRF, **preset_dict)
# Show summary
snl.summary()
[24]:
Type Start End Population ODE Rt theta kappa rho sigma tau alpha1 [-] 1/alpha2 [day] 1/beta [day] 1/gamma [day]
0th Past 02Jan2020 31Jan2020 1000000 SIR-F 2.5 0.002 0.005 0.2 0.075 1440 0.002 200.0 5.0 13.0

Set the 1st phase with the same parameter values.

[25]:
# Add main scenario
snl.add(end_date="31Dec2020", name="Main")
snl.summary()
[25]:
Type Start End Population ODE Rt theta kappa rho sigma tau alpha1 [-] 1/alpha2 [day] 1/beta [day] 1/gamma [day]
0th Past 02Jan2020 31Jan2020 1000000 SIR-F 2.5 0.002 0.005 0.2 0.075 1440 0.002 200.0 5.0 13.0
1st Future 01Feb2020 31Dec2020 1000000 SIR-F 2.5 0.002 0.005 0.2 0.075 1440 0.002 200.0 5.0 13.0

Copy the main scenario and name it as Lockdown scenario. Scenario.clear() removes the future phase (th 1st phase here) and we will register th 1st phase with halved \(\rho\) value. Lockdown is supposed to reduce effective contact rate.

[26]:
# Add lockdown scenario
snl.clear(name="Lockdown")
# Get rho value of the 0th phase and halve it
rho_lock = snl.get("rho", phase="0th") * 0.5
# Add th 1st phase with the calculated rho value
snl.add(end_date="31Dec2020", name="Lockdown", rho=rho_lock)
[26]:
<covsirphy.analysis.scenario.Scenario at 0x7f3a9c9d3ee0>

Next, we define medicine scenario. \(\kappa\) will be halved and \(\sigma\) will be doubled. New medicines may reduce the severity rate and enhance recovery.

[27]:
# Add medicine scenario
snl.clear(name="Medicine")
kappa_med = snl.get("kappa", phase="0th") * 0.5
sigma_med = snl.get("sigma", phase="0th") * 2
snl.add(end_date="31Dec2020", name="Medicine", kappa=kappa_med, sigma=sigma_med)
[27]:
<covsirphy.analysis.scenario.Scenario at 0x7f3a9c9d3ee0>

We define vaccine scenario. As noted in “SIR-F model with vaccination” section, vaccination impacts on \(\sigma\) and \(\kappa\) with depending on its effectivity and safety. If vaccinations impact on infectivity, \(\rho\) value will be also changed.

[28]:
# Add vaccine scenario
snl.clear(name="Vaccine")
rho_vac = snl.get("rho", phase="0th") * 0.8
kappa_vac = snl.get("kappa", phase="0th") * 0.6
sigma_vac = snl.get("sigma", phase="0th") * 1.2
snl.add(end_date="31Dec2020", name="Vaccine",  rho=rho_vac, kappa=kappa_vac, sigma=sigma_vac)
[28]:
<covsirphy.analysis.scenario.Scenario at 0x7f3a9c9d3ee0>

See the phase settings with Scenario.summary() as always.

[29]:
# Show summary
snl.summary()
[29]:
Type Start End Population ODE Rt theta kappa rho sigma tau alpha1 [-] 1/alpha2 [day] 1/beta [day] 1/gamma [day]
Scenario Phase
Main 0th Past 02Jan2020 31Jan2020 1000000 SIR-F 2.50 0.002 0.0050 0.20 0.075 1440 0.002 200.0 5.0 13.0
1st Future 01Feb2020 31Dec2020 1000000 SIR-F 2.50 0.002 0.0050 0.20 0.075 1440 0.002 200.0 5.0 13.0
Lockdown 0th Past 02Jan2020 31Jan2020 1000000 SIR-F 2.50 0.002 0.0050 0.20 0.075 1440 0.002 200.0 5.0 13.0
1st Future 01Feb2020 31Dec2020 1000000 SIR-F 1.25 0.002 0.0050 0.10 0.075 1440 0.002 200.0 10.0 13.0
Medicine 0th Past 02Jan2020 31Jan2020 1000000 SIR-F 2.50 0.002 0.0050 0.20 0.075 1440 0.002 200.0 5.0 13.0
1st Future 01Feb2020 31Dec2020 1000000 SIR-F 1.31 0.002 0.0025 0.20 0.150 1440 0.002 400.0 5.0 6.0
Vaccine 0th Past 02Jan2020 31Jan2020 1000000 SIR-F 2.50 0.002 0.0050 0.20 0.075 1440 0.002 200.0 5.0 13.0
1st Future 01Feb2020 31Dec2020 1000000 SIR-F 1.72 0.002 0.0030 0.16 0.090 1440 0.002 333.0 6.0 11.0

Show the parameter setting with Scenario.history().

[30]:
# Show the history of rho as a dataframe and a figure
# we can set theta/kappa/rho/sigma for SIR-F model
snl.history(target="rho").head()
_images/usage_theoretical_73_0.png
[30]:
Scenario Lockdown Main Medicine Vaccine
Date
2020-01-02 0.2 0.2 0.2 0.2
2020-01-03 0.2 0.2 0.2 0.2
2020-01-04 0.2 0.2 0.2 0.2
2020-01-05 0.2 0.2 0.2 0.2
2020-01-06 0.2 0.2 0.2 0.2

Compare the scenarios

We will compare the scenarios to discuss the impact of changing parameters. We have some methods for that.

  1. Scenario.describe() shows representative values as a dataframe.

  2. Scenario.history(target="Rt") shows the history of automatically calculated reproduction number.

  3. Scenario.history(target="variable name") shows simulated number of cases for the specified variable.

  4. If you have any ideas, please create issues! :)

1. Show representative values

Scenario.describe() compares

  • max number of infected cases and the date with the max value,

  • the number of confirmed/infected/fatal cases on the next date of the last phase, and

  • reproduction numbers for the phases with different numbers.

[31]:
# Describe the scenarios
snl.describe()
[31]:
max(Infected) argmax(Infected) Confirmed on 31Dec2020 Infected on 31Dec2020 Fatal on 31Dec2020 1st_Rt
Main 232800 02Mar2020 890043 0 57295 2.50
Lockdown 44360 12Mar2020 436501 137 28090 1.25
Medicine 51404 24Feb2020 483356 0 9911 1.31
Vaccine 113534 09Mar2020 708848 0 24917 1.72

2. History of reproduction number

Scenario.history(target="Rt") shows the history of reproduction number.

[32]:
# Show the history of reproduction number
_ = snl.history(target="Rt")
_images/usage_theoretical_78_0.png

3. Simulated number of cases of the specified variable

We can also set Confired/Infected/Fatal/Recovered as the target of Scenario.history().

[33]:
# The number of infected cases
_ = snl.history(target="Infected")
_images/usage_theoretical_80_0.png
[34]:
# The number of fatal cases
_ = snl.history(target="Fatal")
_images/usage_theoretical_81_0.png

Simulation of each scenario

We can simulate the all kind of the number of cases for a specific scenario with Scenario.simulate().

[35]:
# Main scenario
_ = snl.simulate(name="Main")
_images/usage_theoretical_83_0.png
[36]:
# Lockdown scenario
_ = snl.simulate(name="Lockdown")
_images/usage_theoretical_84_0.png
[37]:
# Medicine scenario
_ = snl.simulate(name="Medicine")
_images/usage_theoretical_85_0.png
[38]:
# Vaccine scenario
_ = snl.simulate(name="Vaccine")
_images/usage_theoretical_86_0.png

For futher analysis, let’s change the parameter settings and add new scenarios!