Usage: scenario analysis

Open In Colab

This is a quick tour of CovsirPhy. Details scenario analysis will be explained. “Scenario analysis” means that we calculate the number of cases in the future phases with some sets of ODE parameter values. With this analysis, we can estimate the impact of our activities against the outbreak on the number of cases.

Preparation

Prepare the packages.

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

Dataset preparation

Download the datasets to “../input” directory and load them.
Please refer to Usage: datasets for the details.
[2]:
loader = cs.DataLoader("../input")
# The number of cases and population values
jhu_data = loader.jhu()
# Government Response Tracker (OxCGRT)
oxcgrt_data = loader.oxcgrt()
# The number of tests
pcr_data = loader.pcr()
# The number of vaccinations
vaccine_data = loader.vaccine()

Start scenario analysis

As an example, we will analysis the number of cases in Japan. covsirphy.Scenario is the interface for analysis. Please specify the area (country: required, province: optional) when creating the instance and register the datasets with Scenario.register(). As the extra datasets, we can select OxCGRTData, PCRData, VaccineData and CountryData.

[3]:
# Specify country and province (optinal) names
snl = cs.Scenario(country="Japan", province=None)
# Register datasets
snl.register(jhu_data, extras=[oxcgrt_data, pcr_data, vaccine_data])
We call JHUData as “the main datasets” because they are required to calculate the number of susceptible/infected/recovered/fatal cases. These variables are used in SIR-F model.
The other datasets are called as “the extra datasets” and they will be used to predict the futute parameter values of SIR-F model for forecasting the number of cases with some scenarios.
Note:
When we use version <= 2.19.1, instance of PopulationData is required because JHUData does not include population values.
population_data = data_loader.population()
snl.register(jhu_data, population_data)

Additional information:

Display/save figures

We have intactive mode and script mode to display/save figures.

When use use interactive shells, including Jupyter Notebook, we can choose either “interactive shell mode” or “script mode” as follows.

Interactive mode:
Figures will be displayed as the output of code cells.
[4]:
# Choose interactive mode (default is True when we use interactive shells)
snl.interactive = True

When you want to turn-off interactive mode temporarlly, set False as Scenario.interactive or apply show_figure=False as an argument of methods, including Scenario.records(). Methods with figures will be shown later in this tutorial.

[5]:
# apply "show_figures=False" to turn-off intaractive mode temporally
# snl.records(show_figure=False)
Script mode:
In script mode, figures will not be displayed. When filenames were applied to the methods as filename argument, figures will be saved in your local environment.
[6]:
# Stop displaying figures
# snl.interactive = False
# With this mode we can save figures, specifying "filename" argument
# snl.records(filename="records.jpg")

When we run codes as a script (eg. python scenario_analysis.py), only “script mode” is selectable and Scenario.interactive is always False. Figures will be saved when filenames are specified with filename argument.

Because some methods, including Scenario.summary(), return dataframes (pandas.DataFrame), we can save them as CSV files etc. using .to_csv(filename, index=True).

We can produce filenames more easily with Filer class. Please refer to the scripts in example directory of the repository.

[7]:
filer = cs.Filer(directory="output", prefix="jpn", suffix=None, numbering="01")
# filer.png("records")
# -> {"filename": "<absolute path>/output/jpn_01_records.png"}
# filer.jpg("records")
# -> {"filename": "<absolute path>/output/jpn_01_records.jpg"}
# filer.csv("records", index=True)
# -> {"path_or_buf": "<absolute path>/output/jpn_01_records.csv", index: True}

We can save files more easily with Filer as follows.

[8]:
# record_df = snl.records(**filer.png("records"))
# record_df.to_csv(**filer.csv("records", index=False))

Backup/restore scenario

From development version 2.19.1-kappa, we have Scenario.backup(filename) and Scenario.restore(filename) to backup/restore timepoints and phase information. This will be helpful when we perform parameter estimation using a server and simulate with the estimated parameter values using local machines.

Note that we need to execute Scenario.register() in advance to set timepoints with Scenario.restore().

Backup information:

[9]:
# backupfile_dict = cs.Filer(directory="output")
# snl.backup(**backupfile_dict)

Restore information:

[10]:
# backupfile_dict = cs.Filer(directory="output")
# snl = cs.Scenario(country="Japan")
# snl.register(jhu_data)
# snl.restore(**backupfile_dict).summary()

After restoring information, we can skip Scenario.trend() and Scenario.estiate() (their functionarities will be explaiend later).

Check records

Let’s see the records at first. Scenario.records() method return the records as a pandas dataframe and show a line plot. Some kind of complement will be done for analysis, if necessary.

Scenario.records() shows the number of infected/recovered/fatal cases as default. Using variables argument, we can set the variables to show. Here, we chack the number of confirmed/fatal/recovered cases. They are cumurative values.

[11]:
snl.records(variables="CFR").tail()
# This is the same as
# snl.records(variables=["Confirmed", "Fatal", "Recovered"])
_images/usage_quick_30_0.png
[11]:
Date Confirmed Fatal Recovered
487 2021-06-07 762401 12677 706644
488 2021-06-08 763891 12743 710216
489 2021-06-09 765619 12743 713807
490 2021-06-10 767808 13841 717882
491 2021-06-11 769920 13907 721260

The number of infected cases on date:

[12]:
_ = snl.records(variables="I")
# This is the same as
# snl.records(variables=["Infected"])
_images/usage_quick_32_0.png

All available variables can be retrieved with variables="all".

[13]:
df = snl.records(variables="all", show_figure=False)
pprint(df.set_index("Date").columns.tolist(), compact=True)
['Confirmed', 'Infected', 'Fatal', 'Recovered', 'Susceptible', 'School_closing',
 'Workplace_closing', 'Cancel_events', 'Gatherings_restrictions',
 'Transport_closing', 'Stay_home_restrictions',
 'Internal_movement_restrictions', 'International_movement_restrictions',
 'Information_campaigns', 'Testing_policy', 'Contact_tracing',
 'Stringency_index', 'Tests', 'Tests_diff', 'Vaccinations', 'Vaccinated_once',
 'Vaccinated_full']

We can specify the variables to show.

[14]:
snl.records(variables=["Vaccinated_once"]).tail()
_images/usage_quick_36_0.png
[14]:
Date Vaccinated_once
487 2021-06-07 13749628
488 2021-06-08 14505288
489 2021-06-09 15223263
490 2021-06-10 15942264
491 2021-06-11 15942264

We can calculate the number of daily new cases with Scenario.record_diff() method.

[15]:
# Acceptable variables are the same as Scenario.records()
_ = snl.records_diff(variables="C", window=7)
_images/usage_quick_38_0.png

Scenario.show_complement() method is useful to show the kinds of complement. The details of complement are explained in Usage: datasets section.

[16]:
# Show the details of complement
snl.show_complement()
[16]:
Country Province Monotonic_confirmed Monotonic_fatal Monotonic_recovered Full_recovered Partial_recovered
0 Japan - False True True False True

S-R trend analysis

S-R trend analysis finds the change points of SIR-derived ODE parameters. Details will be explained in Usage (details: phases). Phases will be separated with dotted lines. i.e. Dot lines indicate the start dates of phases.

[17]:
snl.trend().summary()
_images/usage_quick_42_0.png
[17]:
Type Start End Population
0th Past 06Feb2020 18Feb2020 126529100
1st Past 19Feb2020 05Mar2020 126529100
2nd Past 06Mar2020 13Mar2020 126529100
3rd Past 14Mar2020 24Mar2020 126529100
4th Past 25Mar2020 02Apr2020 126529100
5th Past 03Apr2020 10Apr2020 126529100
6th Past 11Apr2020 23Apr2020 126529100
7th Past 24Apr2020 02May2020 126529100
8th Past 03May2020 13May2020 126529100
9th Past 14May2020 28May2020 126529100
10th Past 29May2020 10Jun2020 126529100
11th Past 11Jun2020 23Jun2020 126529100
12th Past 24Jun2020 03Jul2020 126529100
13th Past 04Jul2020 16Jul2020 126529100
14th Past 17Jul2020 30Jul2020 126529100
15th Past 31Jul2020 13Aug2020 126529100
16th Past 14Aug2020 26Aug2020 126529100
17th Past 27Aug2020 10Sep2020 126529100
18th Past 11Sep2020 25Sep2020 126529100
19th Past 26Sep2020 09Oct2020 126529100
20th Past 10Oct2020 17Oct2020 126529100
21st Past 18Oct2020 26Oct2020 126529100
22nd Past 27Oct2020 09Nov2020 126529100
23rd Past 10Nov2020 21Nov2020 126529100
24th Past 22Nov2020 03Dec2020 126529100
25th Past 04Dec2020 16Dec2020 126529100
26th Past 17Dec2020 29Dec2020 126529100
27th Past 30Dec2020 11Jan2021 126529100
28th Past 12Jan2021 20Jan2021 126529100
29th Past 21Jan2021 28Jan2021 126529100
30th Past 29Jan2021 05Feb2021 126529100
31st Past 06Feb2021 15Feb2021 126529100
32nd Past 16Feb2021 25Feb2021 126529100
33rd Past 26Feb2021 05Mar2021 126529100
34th Past 06Mar2021 15Mar2021 126529100
35th Past 16Mar2021 23Mar2021 126529100
36th Past 24Mar2021 06Apr2021 126529100
37th Past 07Apr2021 17Apr2021 126529100
38th Past 18Apr2021 30Apr2021 126529100
39th Past 01May2021 14May2021 126529100
40th Past 15May2021 27May2021 126529100
41st Past 28May2021 11Jun2021 126529100

Parameter estimation of ODE models

Here, we will estimate the tau value [min] (using grid search) and parameter values of SIR-derived models using Optuna package (automated hyperparameter optimization framework). As an example, we use SIR-F model. Details of models will be explained in Usage (details: theoritical datasets).

We can select the model from SIR, SIRD and SIR-F model for parameter estimation. SIR-FV model (completely deprecated) and SEWIR-F model cannot be used.

[18]:
# Estimate the tau value and parameter values of SIR-F model
# Default value of timeout in each phase is 180 sec
snl.estimate(cs.SIRF, timeout=180)

<SIR-F model: parameter estimation>
Running optimization with 4 CPUs...
         3rd phase (14Mar2020 - 24Mar2020): finished  863 trials in 0 min 20 sec
         0th phase (06Feb2020 - 18Feb2020): finished 1162 trials in 0 min 25 sec
         6th phase (11Apr2020 - 23Apr2020): finished 1146 trials in 0 min 30 sec
         9th phase (14May2020 - 28May2020): finished 1236 trials in 0 min 35 sec
        10th phase (29May2020 - 10Jun2020): finished  208 trials in 0 min  5 sec
         1st phase (19Feb2020 - 05Mar2020): finished  707 trials in 0 min 20 sec
        11th phase (11Jun2020 - 23Jun2020): finished  438 trials in 0 min 10 sec
         4th phase (25Mar2020 - 02Apr2020): finished 1106 trials in 0 min 35 sec
         7th phase (24Apr2020 - 02May2020): finished  807 trials in 0 min 25 sec
         2nd phase (06Mar2020 - 13Mar2020): finished  931 trials in 0 min 25 sec
        12th phase (24Jun2020 - 03Jul2020): finished 1198 trials in 0 min 35 sec
         8th phase (03May2020 - 13May2020): finished 1251 trials in 0 min 35 sec
         5th phase (03Apr2020 - 10Apr2020): finished 1372 trials in 0 min 40 sec
        15th phase (31Jul2020 - 13Aug2020): finished 1602 trials in 0 min 50 sec
        21st phase (18Oct2020 - 26Oct2020): finished  967 trials in 0 min 25 sec
        16th phase (14Aug2020 - 26Aug2020): finished  466 trials in 0 min 10 sec
        13th phase (04Jul2020 - 16Jul2020): finished 1485 trials in 0 min 45 sec
        18th phase (11Sep2020 - 25Sep2020): finished 1475 trials in 0 min 45 sec
        17th phase (27Aug2020 - 10Sep2020): finished  628 trials in 0 min 20 sec
        19th phase (26Sep2020 - 09Oct2020): finished  902 trials in 0 min 35 sec
        22nd phase (27Oct2020 - 09Nov2020): finished 1519 trials in 1 min  5 sec
        14th phase (17Jul2020 - 30Jul2020): finished 1274 trials in 0 min 55 sec
        23rd phase (10Nov2020 - 21Nov2020): finished  180 trials in 0 min  5 sec
        24th phase (22Nov2020 - 03Dec2020): finished 1116 trials in 0 min 45 sec
        30th phase (29Jan2021 - 05Feb2021): finished  199 trials in 0 min  5 sec
        25th phase (04Dec2020 - 16Dec2020): finished  196 trials in 0 min  5 sec
        31st phase (06Feb2021 - 15Feb2021): finished  181 trials in 0 min  5 sec
        32nd phase (16Feb2021 - 25Feb2021): finished  193 trials in 0 min  5 sec
        33rd phase (26Feb2021 - 05Mar2021): finished  172 trials in 0 min  5 sec
        34th phase (06Mar2021 - 15Mar2021): finished  987 trials in 0 min 40 sec
        20th phase (10Oct2020 - 17Oct2020): finished 1644 trials in 1 min 20 sec
        26th phase (17Dec2020 - 29Dec2020): finished 1256 trials in 0 min 55 sec
        35th phase (16Mar2021 - 23Mar2021): finished  194 trials in 0 min  5 sec
        39th phase (01May2021 - 14May2021): finished  227 trials in 0 min  5 sec
        27th phase (30Dec2020 - 11Jan2021): finished 1628 trials in 1 min 15 sec
        28th phase (12Jan2021 - 20Jan2021): finished  561 trials in 0 min 15 sec
        29th phase (21Jan2021 - 28Jan2021): finished  541 trials in 0 min 15 sec
        36th phase (24Mar2021 - 06Apr2021): finished 1413 trials in 0 min 50 sec
        40th phase (15May2021 - 27May2021): finished 1680 trials in 1 min  0 sec
        37th phase (07Apr2021 - 17Apr2021): finished 1291 trials in 0 min 40 sec
        38th phase (18Apr2021 - 30Apr2021): finished  251 trials in 0 min  5 sec
        41st phase (28May2021 - 11Jun2021): finished 1551 trials in 0 min 50 sec
Completed optimization. Total: 6 min 13 sec
[19]:
# Show the summary of parameter estimation
snl.summary()
[19]:
Type Start End Population ODE Rt theta kappa rho sigma tau alpha1 [-] 1/alpha2 [day] 1/beta [day] 1/gamma [day] RMSLE Trials Runtime
0th Past 06Feb2020 18Feb2020 126529100 SIR-F 2.32 0.027273 0.000000 0.066019 0.027718 1440 - - - - 0.142890 1162.0 0 min 25 sec
1st Past 19Feb2020 05Mar2020 126529100 SIR-F 4.93 0.007714 0.001322 0.136657 0.026198 1440 0.008 756.0 7.0 38.0 0.170717 707.0 0 min 20 sec
2nd Past 06Mar2020 13Mar2020 126529100 SIR-F 5.06 0.022125 0.000656 0.126654 0.023830 1440 0.022 1524.0 7.0 41.0 0.093897 931.0 0 min 25 sec
3rd Past 14Mar2020 24Mar2020 126529100 SIR-F 2.57 0.000096 0.002956 0.063547 0.021736 1440 0.0 338.0 15.0 46.0 0.029797 863.0 0 min 20 sec
4th Past 25Mar2020 02Apr2020 126529100 SIR-F 5.10 0.008478 0.000606 0.117177 0.022166 1440 0.008 1650.0 8.0 45.0 0.036258 1106.0 0 min 35 sec
5th Past 03Apr2020 10Apr2020 126529100 SIR-F 11.67 0.000566 0.001143 0.124003 0.009478 1440 0.001 874.0 8.0 105.0 0.025583 1372.0 0 min 40 sec
6th Past 11Apr2020 23Apr2020 126529100 SIR-F 6.81 0.000500 0.001369 0.076414 0.009854 1440 0.0 730.0 13.0 101.0 0.116075 1146.0 0 min 30 sec
7th Past 24Apr2020 02May2020 126529100 SIR-F 1.44 0.001306 0.001491 0.028691 0.018346 1440 0.001 670.0 34.0 54.0 0.026841 807.0 0 min 25 sec
8th Past 03May2020 13May2020 126529100 SIR-F 0.16 0.048667 0.001391 0.009444 0.054692 1440 0.049 718.0 105.0 18.0 0.093612 1251.0 0 min 35 sec
9th Past 14May2020 28May2020 126529100 SIR-F 0.13 0.114805 0.002110 0.012507 0.085305 1440 0.115 473.0 79.0 11.0 0.022470 1236.0 0 min 35 sec
10th Past 29May2020 10Jun2020 126529100 SIR-F 0.46 0.033851 0.001717 0.032575 0.066187 1440 0.034 582.0 30.0 15.0 0.016150 208.0 0 min 5 sec
11th Past 11Jun2020 23Jun2020 126529100 SIR-F 0.75 0.007468 0.001961 0.063449 0.081788 1440 0.007 509.0 15.0 12.0 0.030209 438.0 0 min 10 sec
12th Past 24Jun2020 03Jul2020 126529100 SIR-F 1.52 0.000097 0.001568 0.109865 0.070674 1440 0.0 637.0 9.0 14.0 0.018046 1198.0 0 min 35 sec
13th Past 04Jul2020 16Jul2020 126529100 SIR-F 2.10 0.000240 0.000338 0.138912 0.065806 1440 0.0 2958.0 7.0 15.0 0.022398 1485.0 0 min 45 sec
14th Past 17Jul2020 30Jul2020 126529100 SIR-F 2.00 0.000271 0.000185 0.135558 0.067589 1440 0.0 5400.0 7.0 14.0 0.021450 1274.0 0 min 55 sec
15th Past 31Jul2020 13Aug2020 126529100 SIR-F 1.67 0.000058 0.000327 0.122194 0.072657 1440 0.0 3056.0 8.0 13.0 0.046245 1602.0 0 min 50 sec
16th Past 14Aug2020 26Aug2020 126529100 SIR-F 0.83 0.000252 0.000794 0.079969 0.094974 1440 0.0 1259.0 12.0 10.0 0.012330 466.0 0 min 10 sec
17th Past 27Aug2020 10Sep2020 126529100 SIR-F 0.72 0.000429 0.001279 0.071445 0.097862 1440 0.0 781.0 13.0 10.0 0.010068 628.0 0 min 20 sec
18th Past 11Sep2020 25Sep2020 126529100 SIR-F 0.87 0.009109 0.000541 0.077862 0.088208 1440 0.009 1849.0 12.0 11.0 0.017670 1475.0 0 min 45 sec
19th Past 26Sep2020 09Oct2020 126529100 SIR-F 1.00 0.000173 0.001029 0.088567 0.087578 1440 0.0 972.0 11.0 11.0 0.012229 902.0 0 min 35 sec
20th Past 10Oct2020 17Oct2020 126529100 SIR-F 1.04 0.000024 0.000780 0.092602 0.087889 1440 0.0 1281.0 10.0 11.0 0.014964 1644.0 1 min 20 sec
21st Past 18Oct2020 26Oct2020 126529100 SIR-F 0.97 0.000175 0.000894 0.100391 0.102726 1440 0.0 1118.0 9.0 9.0 0.020754 967.0 0 min 25 sec
22nd Past 27Oct2020 09Nov2020 126529100 SIR-F 1.32 0.000480 0.001075 0.113700 0.085210 1440 0.0 929.0 8.0 11.0 0.016758 1519.0 1 min 5 sec
23rd Past 10Nov2020 21Nov2020 126529100 SIR-F 1.58 0.001401 0.000735 0.142531 0.089508 1440 0.001 1360.0 7.0 11.0 0.010060 180.0 0 min 5 sec
24th Past 22Nov2020 03Dec2020 126529100 SIR-F 1.29 0.000446 0.000861 0.105665 0.081151 1440 0.0 1161.0 9.0 12.0 0.014784 1116.0 0 min 45 sec
25th Past 04Dec2020 16Dec2020 126529100 SIR-F 1.18 0.000290 0.001435 0.096849 0.080680 1440 0.0 696.0 10.0 12.0 0.011146 196.0 0 min 5 sec
26th Past 17Dec2020 29Dec2020 126529100 SIR-F 1.20 0.000062 0.001494 0.102799 0.083850 1440 0.0 669.0 9.0 11.0 0.011793 1256.0 0 min 55 sec
27th Past 30Dec2020 11Jan2021 126529100 SIR-F 1.62 0.000261 0.001226 0.106819 0.064709 1440 0.0 815.0 9.0 15.0 0.026075 1628.0 1 min 15 sec
28th Past 12Jan2021 20Jan2021 126529100 SIR-F 1.25 0.000405 0.000916 0.089477 0.070369 1440 0.0 1091.0 11.0 14.0 0.008139 561.0 0 min 15 sec
29th Past 21Jan2021 28Jan2021 126529100 SIR-F 0.77 0.000124 0.001285 0.069112 0.088504 1440 0.0 778.0 14.0 11.0 0.014629 541.0 0 min 15 sec
30th Past 29Jan2021 05Feb2021 126529100 SIR-F 0.61 0.000417 0.002017 0.064249 0.102909 1440 0.0 495.0 15.0 9.0 0.015993 199.0 0 min 5 sec
31st Past 06Feb2021 15Feb2021 126529100 SIR-F 0.53 0.009714 0.001959 0.054271 0.099168 1440 0.01 510.0 18.0 10.0 0.007579 181.0 0 min 5 sec
32nd Past 16Feb2021 25Feb2021 126529100 SIR-F 0.60 0.001622 0.003221 0.065833 0.105507 1440 0.002 310.0 15.0 9.0 0.020170 193.0 0 min 5 sec
33rd Past 26Feb2021 05Mar2021 126529100 SIR-F 0.72 0.000158 0.003500 0.067449 0.090271 1440 0.0 285.0 14.0 11.0 0.006750 172.0 0 min 5 sec
34th Past 06Mar2021 15Mar2021 126529100 SIR-F 1.00 0.013967 0.002137 0.087264 0.083969 1440 0.014 467.0 11.0 11.0 0.010843 987.0 0 min 40 sec
35th Past 16Mar2021 23Mar2021 126529100 SIR-F 1.16 0.009626 0.001781 0.102443 0.085748 1440 0.01 561.0 9.0 11.0 0.009516 194.0 0 min 5 sec
36th Past 24Mar2021 06Apr2021 126529100 SIR-F 1.63 0.000315 0.001435 0.114177 0.068588 1440 0.0 697.0 8.0 14.0 0.010445 1413.0 0 min 50 sec
37th Past 07Apr2021 17Apr2021 126529100 SIR-F 1.52 0.000556 0.000809 0.114059 0.074379 1440 0.001 1235.0 8.0 13.0 0.016669 1291.0 0 min 40 sec
38th Past 18Apr2021 30Apr2021 126529100 SIR-F 1.50 0.001002 0.000776 0.096300 0.063208 1440 0.001 1288.0 10.0 15.0 0.008988 251.0 0 min 5 sec
39th Past 01May2021 14May2021 126529100 SIR-F 1.29 0.000116 0.001031 0.085381 0.065184 1440 0.0 970.0 11.0 15.0 0.007359 227.0 0 min 5 sec
40th Past 15May2021 27May2021 126529100 SIR-F 0.83 0.000110 0.001374 0.071418 0.084350 1440 0.0 727.0 14.0 11.0 0.010721 1680.0 1 min 0 sec
41st Past 28May2021 11Jun2021 126529100 SIR-F 0.59 0.000357 0.001973 0.055327 0.092367 1440 0.0 506.0 18.0 10.0 0.016189 1551.0 0 min 50 sec

Evaluation of estimation accuracy

Accuracy of parameter estimation can be evaluated with RMSLE (Root Mean Squared Log Error) score.

\begin{align*} \mathrm{RMSLE} = \sqrt{\cfrac{1}{n}\sum_{i=1}^{n}(log_{10}(A_{i} + 1) - log_{10}(P_{i} + 1))^2} \end{align*}

Where \(A\) is the observed (actual) values, \(P\) is estimated (predicted) values. Variables are \(S (i=1), I (i=2), R (i=3)\ \mathrm{and}\ F (i=n=4)\) for SIR-F model. When RMSLE score is low, hyperparameter estimation is highly accurate. Please refer to external sites, including Medium: What’s the Difference Between RMSE and RMSLE?

[20]:
# Show RMSLE scores with the number of optimization trials and runtime for phases
snl.summary(columns=["Start", "End", "RMSLE", "Trials", "Runtime"])
[20]:
Start End RMSLE Trials Runtime
0th 06Feb2020 18Feb2020 0.142890 1162.0 0 min 25 sec
1st 19Feb2020 05Mar2020 0.170717 707.0 0 min 20 sec
2nd 06Mar2020 13Mar2020 0.093897 931.0 0 min 25 sec
3rd 14Mar2020 24Mar2020 0.029797 863.0 0 min 20 sec
4th 25Mar2020 02Apr2020 0.036258 1106.0 0 min 35 sec
5th 03Apr2020 10Apr2020 0.025583 1372.0 0 min 40 sec
6th 11Apr2020 23Apr2020 0.116075 1146.0 0 min 30 sec
7th 24Apr2020 02May2020 0.026841 807.0 0 min 25 sec
8th 03May2020 13May2020 0.093612 1251.0 0 min 35 sec
9th 14May2020 28May2020 0.022470 1236.0 0 min 35 sec
10th 29May2020 10Jun2020 0.016150 208.0 0 min 5 sec
11th 11Jun2020 23Jun2020 0.030209 438.0 0 min 10 sec
12th 24Jun2020 03Jul2020 0.018046 1198.0 0 min 35 sec
13th 04Jul2020 16Jul2020 0.022398 1485.0 0 min 45 sec
14th 17Jul2020 30Jul2020 0.021450 1274.0 0 min 55 sec
15th 31Jul2020 13Aug2020 0.046245 1602.0 0 min 50 sec
16th 14Aug2020 26Aug2020 0.012330 466.0 0 min 10 sec
17th 27Aug2020 10Sep2020 0.010068 628.0 0 min 20 sec
18th 11Sep2020 25Sep2020 0.017670 1475.0 0 min 45 sec
19th 26Sep2020 09Oct2020 0.012229 902.0 0 min 35 sec
20th 10Oct2020 17Oct2020 0.014964 1644.0 1 min 20 sec
21st 18Oct2020 26Oct2020 0.020754 967.0 0 min 25 sec
22nd 27Oct2020 09Nov2020 0.016758 1519.0 1 min 5 sec
23rd 10Nov2020 21Nov2020 0.010060 180.0 0 min 5 sec
24th 22Nov2020 03Dec2020 0.014784 1116.0 0 min 45 sec
25th 04Dec2020 16Dec2020 0.011146 196.0 0 min 5 sec
26th 17Dec2020 29Dec2020 0.011793 1256.0 0 min 55 sec
27th 30Dec2020 11Jan2021 0.026075 1628.0 1 min 15 sec
28th 12Jan2021 20Jan2021 0.008139 561.0 0 min 15 sec
29th 21Jan2021 28Jan2021 0.014629 541.0 0 min 15 sec
30th 29Jan2021 05Feb2021 0.015993 199.0 0 min 5 sec
31st 06Feb2021 15Feb2021 0.007579 181.0 0 min 5 sec
32nd 16Feb2021 25Feb2021 0.020170 193.0 0 min 5 sec
33rd 26Feb2021 05Mar2021 0.006750 172.0 0 min 5 sec
34th 06Mar2021 15Mar2021 0.010843 987.0 0 min 40 sec
35th 16Mar2021 23Mar2021 0.009516 194.0 0 min 5 sec
36th 24Mar2021 06Apr2021 0.010445 1413.0 0 min 50 sec
37th 07Apr2021 17Apr2021 0.016669 1291.0 0 min 40 sec
38th 18Apr2021 30Apr2021 0.008988 251.0 0 min 5 sec
39th 01May2021 14May2021 0.007359 227.0 0 min 5 sec
40th 15May2021 27May2021 0.010721 1680.0 1 min 0 sec
41st 28May2021 11Jun2021 0.016189 1551.0 0 min 50 sec

Additionally, we can visualize the accuracy with Scenario.estimate_accuracy(), specifing phase name.

[21]:
# Visualize the accuracy for the 2nd phase
snl.estimate_accuracy(phase="2nd")
# phase="last" means the last phases
# snl.estimate_accuracy(phase="last")
_images/usage_quick_49_0.png

We can calculate total score for all phases using Scenario.score() method. Metrics can be selected from MAE, MSE, MSLE, RMSE and RMSLE.

[22]:
# Get total score: metrics="MAE", "MSE", "MSLE", "RMSE" or "RMSLE"
# snl.score(metrics="RMSLE")
metrics_list = ["MAE", "MSE", "MSLE", "RMSE", "RMSLE"]
for metrics in metrics_list:
    metrics_name = metrics.rjust(len(max(metrics_list, key=len)))
    print(f"{metrics_name}: {snl.score(metrics=metrics):.3f}")
  MAE: 731.091
  MSE: 2162748.018
 MSLE: 0.003
 RMSE: 1150.362
RMSLE: 0.058

Get parameter value

We can get the parameter values of a phase using Scenario.get() method.

[23]:
# Get parameter values
snl.get("Rt", phase="4th")
[23]:
5.1
[24]:
# phase="last" means the last phases
snl.get("Rt", phase="last")
[24]:
0.59

Show parameter history

We can get the history of parameter values with a dataframe and a figure.

[25]:
# Get the parameter values as a dataframe
snl.summary(columns=[*cs.SIRF.PARAMETERS, "Rt"])
[25]:
theta kappa rho sigma Rt
0th 0.027273 0.000000 0.066019 0.027718 2.32
1st 0.007714 0.001322 0.136657 0.026198 4.93
2nd 0.022125 0.000656 0.126654 0.023830 5.06
3rd 0.000096 0.002956 0.063547 0.021736 2.57
4th 0.008478 0.000606 0.117177 0.022166 5.10
5th 0.000566 0.001143 0.124003 0.009478 11.67
6th 0.000500 0.001369 0.076414 0.009854 6.81
7th 0.001306 0.001491 0.028691 0.018346 1.44
8th 0.048667 0.001391 0.009444 0.054692 0.16
9th 0.114805 0.002110 0.012507 0.085305 0.13
10th 0.033851 0.001717 0.032575 0.066187 0.46
11th 0.007468 0.001961 0.063449 0.081788 0.75
12th 0.000097 0.001568 0.109865 0.070674 1.52
13th 0.000240 0.000338 0.138912 0.065806 2.10
14th 0.000271 0.000185 0.135558 0.067589 2.00
15th 0.000058 0.000327 0.122194 0.072657 1.67
16th 0.000252 0.000794 0.079969 0.094974 0.83
17th 0.000429 0.001279 0.071445 0.097862 0.72
18th 0.009109 0.000541 0.077862 0.088208 0.87
19th 0.000173 0.001029 0.088567 0.087578 1.00
20th 0.000024 0.000780 0.092602 0.087889 1.04
21st 0.000175 0.000894 0.100391 0.102726 0.97
22nd 0.000480 0.001075 0.113700 0.085210 1.32
23rd 0.001401 0.000735 0.142531 0.089508 1.58
24th 0.000446 0.000861 0.105665 0.081151 1.29
25th 0.000290 0.001435 0.096849 0.080680 1.18
26th 0.000062 0.001494 0.102799 0.083850 1.20
27th 0.000261 0.001226 0.106819 0.064709 1.62
28th 0.000405 0.000916 0.089477 0.070369 1.25
29th 0.000124 0.001285 0.069112 0.088504 0.77
30th 0.000417 0.002017 0.064249 0.102909 0.61
31st 0.009714 0.001959 0.054271 0.099168 0.53
32nd 0.001622 0.003221 0.065833 0.105507 0.60
33rd 0.000158 0.003500 0.067449 0.090271 0.72
34th 0.013967 0.002137 0.087264 0.083969 1.00
35th 0.009626 0.001781 0.102443 0.085748 1.16
36th 0.000315 0.001435 0.114177 0.068588 1.63
37th 0.000556 0.000809 0.114059 0.074379 1.52
38th 0.001002 0.000776 0.096300 0.063208 1.50
39th 0.000116 0.001031 0.085381 0.065184 1.29
40th 0.000110 0.001374 0.071418 0.084350 0.83
41st 0.000357 0.001973 0.055327 0.092367 0.59

Scenario.history() method shows the trajectories of parameters (and the number of cases).

[26]:
_ = snl.history(target="theta", show_legend=False)
_images/usage_quick_58_0.png
[27]:
_ = snl.history(target="kappa", show_legend=False)
_images/usage_quick_59_0.png
[28]:
_ = snl.history(target="rho", show_legend=False)
_images/usage_quick_60_0.png
[29]:
_ = snl.history(target="sigma", show_legend=False)
_images/usage_quick_61_0.png
Notes on the history of \(\sigma\) value in japan (last updated: 28Dec2020):
In Japan, we experienced two waves and we are in third wave. In the first wave (Apr - May), recovery period was too long because collapse of the medical care system occurred and no medicines were found.

Sigma values: the first wave < the second wave > the third wave

However, in the second wave (Jul - Oct), recovery period appears short because we have some effective medicines (not approved, in clinical study), yonger people (people un-associated to sever diseases) were infected.

In the third wave (Nov - ), older people tend to be infected and we are facing with medical collapse at this time…

Show the history of reproduction number

\(R_0\) (“R naught”) means “the average number of secondary infections caused by an infected host” (Infection Modeling — Part 1). When this value is larger than 1, the infection disease is outbreaking.

[30]:
_ = snl.history(target="Rt", show_legend=False)
_images/usage_quick_64_0.png

Simulate the number of cases

We can compare the actual and simulated (with estimated parameter values) number of confirmed/infected/recovered/fatal cases using Scenario.history() method.

[31]:
# Compare the actual values and the main scenario
_ = snl.history("Infected")
_images/usage_quick_66_0.png

When we want to show only one scenario with all variables, we use Scenario.simulate() method.

[32]:
_ =snl.simulate(name="Main")
_images/usage_quick_68_0.png

Main scenario

To investigate the effect of parameter changes, we will perform scenario analysis. In the main scenario, we will assume that the parameter values do not change after the last past phase.

i.e. If the parameter velues will not be changed until 31Jul2021, how many cases will be? We call this scenario as “Main” scenario.

[33]:
# Clear future phases in Main scenario
snl.clear(name="Main")
# Add one future phase 30 days with the parameter set of the last past phase
snl.add(days=30, name="Main")
# Add one future phase until 31Jul2021 with the same parameter set
snl.add(end_date="31Jul2021", name="Main")
# Simulate the number of cases
snl.simulate(name="Main").tail()
_images/usage_quick_70_0.png
[33]:
Date Infected Fatal Recovered
524 2021-07-27 5641 14741 792614
525 2021-07-28 5423 14752 793125
526 2021-07-29 5214 14762 793615
527 2021-07-30 5013 14772 794087
528 2021-07-31 4819 14782 794541

Medicine scenario

To investigate the effect of new medicines, we will assume that \(\sigma\) will be changed in the future phases.

If \(\sigma\) will be 1.2 times in 30 days, how many cases will be? We will call this scenario as “Medicine” scenario.

[34]:
# Calcuate the current sigma value of the last phase
sigma_current = snl.get("sigma", name="Main", phase="last")
sigma_current
[34]:
0.09236710845320784
[35]:
# Sigma value will be double
sigma_new = sigma_current * 1.2
sigma_new
[35]:
0.1108405301438494
[36]:
# Initialize "Medicine" scenario (with the same past phases as that of Main scenario)
snl.clear(name="Medicine")
# Add 30 days as a new future phases with the same parameter set
snl.add(name="Medicine", days=30, sigma=sigma_current)
# Add a phase with doubled sigma value and the same end date with main date
snl.add(name="Medicine", end_date="31Jul2021", sigma=sigma_new)
[36]:
<covsirphy.analysis.scenario.Scenario at 0x7f52955f52e0>

Check summary of future phases.

[37]:
df = snl.summary()
df.loc[df["Type"] == "Future"]
[37]:
Type Start End Population ODE Rt theta kappa rho sigma tau alpha1 [-] 1/alpha2 [day] 1/beta [day] 1/gamma [day] RMSLE Trials Runtime
Scenario Phase
Main 42nd Future 12Jun2021 11Jul2021 126529100 SIR-F 0.59 0.000357 0.001973 0.055327 0.092367 1440 0.0 506.0 18.0 10.0 - - -
43rd Future 12Jul2021 31Jul2021 126529100 SIR-F 0.59 0.000357 0.001973 0.055327 0.092367 1440 0.0 506.0 18.0 10.0 - - -
Medicine 42nd Future 12Jun2021 11Jul2021 126529100 SIR-F 0.59 0.000357 0.001973 0.055327 0.092367 1440 0.0 506.0 18.0 10.0 - - -
43rd Future 12Jul2021 31Jul2021 126529100 SIR-F 0.49 0.000357 0.001973 0.055327 0.110841 1440 0.0 506.0 18.0 9.0 - - -

Simulate the number of cases.

[38]:
_ = snl.simulate(name="Medicine").tail()
_images/usage_quick_78_0.png

Short-term prediction of parameter values

With extra datasets, we can predict the ODE parameter values in the future phases because OxCGRT indicators (policy measures), vaccinations and so on impact on parameter values with delay period (the number of days that indicators take for the ODE parameter values to be changed).

OxCGRT indicators are

  • school_closing,

  • workplace_closing,

  • cancel_events,

  • gatherings_restrictions,

  • transport_closing,

  • stay_home_restrictions,

  • internal_movement_restrictions,

  • international_movement_restrictions,

  • information_campaigns,

  • testing_policy, and

  • contact_tracing.

Scenario.fit() method learns the relationship of indicators (X) and the parameter values (y) with regresion models. X was registered with Scenario.register() and y was calculated with Scenario.estimate() in advance respectively.

Regression models:

  • Indicators -> Parameters with Elastic Net

  • Indicators -> Parameters with Decision Tree Regressor

  • Indicators(n)/Indicators(n-1) -> Parameters(n)/Parameters(n-1) with Elastic Net

  • Indicators(n)/Indicators(n-1) -> Parameters(n)/Parameters(n-1) with Decision Tree Regressor

Scenario.fit() has delay argument and we can apply the following values.

  • delay=20: Exact value of delay period [days].

  • delay=(7, 31): Delay period will be found by grid search in the value range of (7, 31) to minimize test score of learing. This is the default value from version 2.20.0.

  • delay=None: Delay period will be calculated with Scenario.estimate_delay() method. In this method, change point analysis will be done with one selected indicator (“Stringency_index” as default) and variable (“Confirmed” as default) to estimate delay period. This is the default value to version 2.19.1.

When we use delay=None, the indicator used to estimate delay period (“Stringency_index” as default) will be automatically removed from X dataset. If you have some indicators to be removed additionally, please use removed_cols argument (list of indicator names) of Scenario.fit().

As an example, we only use he number of vaccinated (once/full) people and OxCGRT indicators (except for “Stringency_index” because this is calculated with the other OxCGRT indicators).

[39]:
# Create Forecast scenario (copy Main scenario and delete future phases)
snl.clear(name="Forecast", template="Main")
# Fitting with linear regression model (Elastic Net regression)
fit_dict = snl.fit(
    name="Forecast", metric="R2",
    removed_cols=["Stringency_index", "Tests", "Tests_diff", "Vaccinations"]
)

Show determination coeefients of training/test dataset.

[40]:
print(f"R-squared: {fit_dict['score_train']:.3f} (train)")
print(f"R-squared: {fit_dict['score_test']:.3f} (test)")
R-squared: -2.422 (train)
R-squared: 0.000 (test)

Scenario.predict() predicts the parameter values of future phases.

[41]:
# Short-term prediction
snl.predict(name="Forecast")
# We can select list of days to predict optionally
# snl.predict(days=[1, 4], name="Forecast")
[41]:
<covsirphy.analysis.scenario.Scenario at 0x7f52955f52e0>

Show summary of future phases.

[42]:
df = snl.summary()
df.loc[df["Type"] == "Future"]
[42]:
Type Start End Population ODE Rt theta kappa rho sigma tau alpha1 [-] 1/alpha2 [day] 1/beta [day] 1/gamma [day] RMSLE Trials Runtime
Scenario Phase
Main 42nd Future 12Jun2021 11Jul2021 126529100 SIR-F 0.59 0.000357 0.001973 0.055327 0.092367 1440 0.0 506.0 18.0 10.0 - - -
43rd Future 12Jul2021 31Jul2021 126529100 SIR-F 0.59 0.000357 0.001973 0.055327 0.092367 1440 0.0 506.0 18.0 10.0 - - -
Medicine 42nd Future 12Jun2021 11Jul2021 126529100 SIR-F 0.59 0.000357 0.001973 0.055327 0.092367 1440 0.0 506.0 18.0 10.0 - - -
43rd Future 12Jul2021 31Jul2021 126529100 SIR-F 0.49 0.000357 0.001973 0.055327 0.110841 1440 0.0 506.0 18.0 9.0 - - -
Forecast 42nd Future 12Jun2021 12Jul2021 126529100 SIR-F 0.59 0.000357 0.001973 0.055330 0.092370 1440 0.0 506.0 18.0 10.0 - - -

Or, short-cut .fit() and .predict() with Scenario.fit_predict(). We can apply all arguments (of .fit() and .predict()) to Scenario.fit_predict().

[43]:
# Or, when you do not need 'fit_dict',
# snl.fit_predict(name="Forecast").summary(name="Forecast")

To compare this scenario with the other scenarios, we should adjust the last end date with Scenario.adjust_end() because the last end date is differenct from the other scenarios at this time.

[44]:
# Adjust the last end dates
snl.adjust_end()
# Show the last phases of all scenarios
all_df = snl.summary().reset_index()
for name in all_df["Scenario"].unique():
    df = snl.summary(name=name)
    last_end_date = df.loc[df.index[-1], "End"]
    print(f"{name} scenario: to {last_end_date}")
Main scenario: to 31Jul2021
Medicine scenario: to 31Jul2021
Forecast scenario: to 31Jul2021

Simulate the number of cases of forecast scenario.

[45]:
_ = snl.simulate(variables="CFR", name="Forecast").tail()
_images/usage_quick_95_0.png

Compare the number of cases of the all scenario with Scenario.history() and variable name.

[46]:
_ = snl.history("Infected")
_images/usage_quick_97_0.png

We can focus on the values in specified date range with the following arguments.

  • dates: tuple of start date and end date

  • past_days (integer): how many past days to use in calculation from today (Scenario.today property)

  • phases (list of str): phase names to use in calculation

These arguments are effective with Scenario.history(), Scenario.simulate(), Scenario.track() and Scenario.score().

[47]:
# Get the minimum value (from today to future) to set lower limit of y-axis
lower_limit = snl.history("Infected", dates=(snl.today, None), show_figure=False).min().min()
# From today to future (no limitation regarding end date)
_ = snl.history("Infected", dates=(snl.today, None), ylim=(lower_limit, None))
_images/usage_quick_99_0.png

In the past 20 days. Reference date is today (Scenario.today property).

[48]:
_ = snl.history("Infected", past_days=20)
_images/usage_quick_101_0.png

In the selected phases. Here, we will show the 3rd, 4th and 5th phase.

[49]:
_ = snl.history("Infected", phases=["3rd", "4th", "5th"])
_images/usage_quick_103_0.png

Compare the scenarios

We will compare the scenarios with representative values, reproduction number and parameter values. Currently, we can compare the scenarios with the following indexes.

  • max(Infected): max value of Infected

  • argmax(Infected): the date when Infected shows max value

  • Infected on …: Infected on the end date of the last phase

  • Fatal on …: Fatal on the end date of the last phase

[50]:
snl.describe()
[50]:
max(Infected) argmax(Infected) Confirmed on 31Jul2021 Infected on 31Jul2021 Fatal on 31Jul2021 43rd_Rt
Main 76460 14May2021 814142 4819 14782 0.59
Medicine 76460 14May2021 812979 3335 14740 0.49
Forecast 76460 14May2021 814145 4819 14782 0.59
[51]:
_ = snl.history(target="Infected")
_images/usage_quick_106_0.png
[52]:
_ = snl.history(target="Rt")
_images/usage_quick_107_0.png
[53]:
_ = snl.history(target="rho")
_images/usage_quick_108_0.png
[54]:
_ = snl.history(target="sigma")
_images/usage_quick_109_0.png
[55]:
_ = snl.history(target="theta")
_images/usage_quick_110_0.png
[56]:
_ = snl.history(target="kappa")
_images/usage_quick_111_0.png

Change rate of parameters in main scenario

History of each parameter will be shown. Values will be divided by the values in 0th phase.

[57]:
_ = snl.history_rate(name="Main")
_images/usage_quick_113_0.png

Retrospective analysis

We can evaluate the impact of measures using past records. How many people were infected if the parameter values have not changed sinse 01Sep2020?

[58]:
# Perform retrospective analysis
snl_retro = cs.Scenario(country="Japan")
snl_retro.register(jhu_data)
snl_retro.retrospective(
    "01Jan2021", model=cs.SIRF, control="Main", target="Retrospective", timeout=10)

<SIR-F model: parameter estimation>
Running optimization with 4 CPUs...
         9th phase (14May2020 - 28May2020): finished  146 trials in 0 min  5 sec
         3rd phase (14Mar2020 - 24Mar2020): finished  272 trials in 0 min 10 sec
         0th phase (06Feb2020 - 18Feb2020): finished  331 trials in 0 min 10 sec
         6th phase (11Apr2020 - 23Apr2020): finished  266 trials in 0 min 10 sec
        10th phase (29May2020 - 10Jun2020): finished  144 trials in 0 min  5 sec
         1st phase (19Feb2020 - 05Mar2020): finished  278 trials in 0 min 10 sec
         7th phase (24Apr2020 - 02May2020): finished  290 trials in 0 min 10 sec
        11th phase (11Jun2020 - 23Jun2020): finished  282 trials in 0 min 10 sec
         4th phase (25Mar2020 - 02Apr2020): finished  276 trials in 0 min 10 sec
         8th phase (03May2020 - 13May2020): finished  301 trials in 0 min 10 sec         2nd phase (06Mar2020 - 13Mar2020): finished  270 trials in 0 min 10 sec

        12th phase (24Jun2020 - 03Jul2020): finished  271 trials in 0 min 10 sec         5th phase (03Apr2020 - 10Apr2020): finished  301 trials in 0 min 10 sec

        15th phase (31Jul2020 - 13Aug2020): finished  287 trials in 0 min 10 sec
        18th phase (11Sep2020 - 25Sep2020): finished  303 trials in 0 min 10 sec
        13th phase (04Jul2020 - 16Jul2020): finished  296 trials in 0 min 10 sec        21st phase (18Oct2020 - 26Oct2020): finished  290 trials in 0 min 10 sec

        16th phase (14Aug2020 - 26Aug2020): finished  275 trials in 0 min 10 sec
        19th phase (26Sep2020 - 09Oct2020): finished  270 trials in 0 min 10 sec        14th phase (17Jul2020 - 30Jul2020): finished  262 trials in 0 min 10 sec

        22nd phase (27Oct2020 - 09Nov2020): finished  262 trials in 0 min 10 sec
        23rd phase (10Nov2020 - 21Nov2020): finished  140 trials in 0 min  5 sec
        17th phase (27Aug2020 - 10Sep2020): finished  277 trials in 0 min 10 sec
        24th phase (22Nov2020 - 03Dec2020): finished  280 trials in 0 min 10 sec
        20th phase (10Oct2020 - 17Oct2020): finished  273 trials in 0 min 10 sec
        27th phase (30Dec2020 - 31Dec2020): finished  284 trials in 0 min  5 sec
        33rd phase (16Feb2021 - 25Feb2021): finished  145 trials in 0 min  5 sec
        30th phase (21Jan2021 - 28Jan2021): finished  261 trials in 0 min 10 sec
        25th phase (04Dec2020 - 16Dec2020): finished  267 trials in 0 min 10 sec
        34th phase (26Feb2021 - 05Mar2021): finished  142 trials in 0 min  5 sec
        28th phase (01Jan2021 - 11Jan2021): finished  267 trials in 0 min 10 sec
        31st phase (29Jan2021 - 05Feb2021): finished  149 trials in 0 min  5 sec
        32nd phase (06Feb2021 - 15Feb2021): finished  161 trials in 0 min  5 sec
        26th phase (17Dec2020 - 29Dec2020): finished  285 trials in 0 min 10 sec
        35th phase (06Mar2021 - 15Mar2021): finished  285 trials in 0 min 10 sec
        29th phase (12Jan2021 - 20Jan2021): finished  290 trials in 0 min 10 sec
        36th phase (16Mar2021 - 23Mar2021): finished  154 trials in 0 min  5 sec
        39th phase (18Apr2021 - 30Apr2021): finished  151 trials in 0 min  5 sec
        40th phase (01May2021 - 14May2021): finished  171 trials in 0 min  5 sec
        42nd phase (28May2021 - 11Jun2021): finished  308 trials in 0 min 10 sec
        37th phase (24Mar2021 - 06Apr2021): finished  354 trials in 0 min 10 sec
        41st phase (15May2021 - 27May2021): finished  362 trials in 0 min 10 sec
        38th phase (07Apr2021 - 17Apr2021): finished  365 trials in 0 min 10 sec
Completed optimization. Total: 1 min 46 sec

<SIR-F model: parameter estimation>
Running optimization with 4 CPUs...
         2nd phase (06Mar2020 - 13Mar2020): finished  246 trials in 0 min 10 sec
         0th phase (06Feb2020 - 18Feb2020): finished  280 trials in 0 min 10 sec
         6th phase (11Apr2020 - 23Apr2020): finished  244 trials in 0 min 10 sec
         4th phase (25Mar2020 - 02Apr2020): finished  244 trials in 0 min 10 sec
         3rd phase (14Mar2020 - 24Mar2020): finished  260 trials in 0 min 10 sec
         1st phase (19Feb2020 - 05Mar2020): finished  259 trials in 0 min 10 sec         7th phase (24Apr2020 - 02May2020): finished  261 trials in 0 min 10 sec

         5th phase (03Apr2020 - 10Apr2020): finished  259 trials in 0 min 10 sec
        10th phase (29May2020 - 10Jun2020): finished  137 trials in 0 min  5 sec
         8th phase (03May2020 - 13May2020): finished  267 trials in 0 min 10 sec
        12th phase (24Jun2020 - 03Jul2020): finished  255 trials in 0 min 10 sec
        14th phase (17Jul2020 - 30Jul2020): finished  253 trials in 0 min 10 sec
         9th phase (14May2020 - 28May2020): finished  127 trials in 0 min  5 sec
        11th phase (11Jun2020 - 23Jun2020): finished  252 trials in 0 min 10 sec
        15th phase (31Jul2020 - 13Aug2020): finished  234 trials in 0 min 10 sec
        13th phase (04Jul2020 - 16Jul2020): finished  240 trials in 0 min 10 sec
        16th phase (14Aug2020 - 26Aug2020): finished  240 trials in 0 min 10 sec
        18th phase (11Sep2020 - 25Sep2020): finished  247 trials in 0 min 10 sec
        20th phase (10Oct2020 - 17Oct2020): finished  257 trials in 0 min 10 sec
        22nd phase (27Oct2020 - 09Nov2020): finished  262 trials in 0 min 10 sec
        17th phase (27Aug2020 - 10Sep2020): finished  232 trials in 0 min 10 sec
        19th phase (26Sep2020 - 09Oct2020): finished  240 trials in 0 min 10 sec
        23rd phase (10Nov2020 - 21Nov2020): finished  116 trials in 0 min  5 sec
        21st phase (18Oct2020 - 26Oct2020): finished  252 trials in 0 min 10 sec
        24th phase (22Nov2020 - 03Dec2020): finished  306 trials in 0 min 10 sec
        26th phase (17Dec2020 - 29Dec2020): finished  309 trials in 0 min 10 sec
        28th phase (01Jan2021 - 11Jun2021): finished  292 trials in 0 min 10 sec
        27th phase (30Dec2020 - 31Dec2020): finished  312 trials in 0 min  5 sec
        25th phase (04Dec2020 - 16Dec2020): finished  377 trials in 0 min 10 sec
Completed optimization. Total: 1 min 16 sec
[59]:
# Show the summary of estimation
cols = ["Start", "End", "ODE", "Rt", *cs.SIRF.PARAMETERS] + ["RMSLE", "Trials", "Runtime"]
snl_retro.summary(columns=cols)
[59]:
Start End ODE Rt theta kappa rho sigma RMSLE Trials Runtime
Scenario Phase
Main 0th 06Feb2020 18Feb2020 SIR-F 2.23 0.027740 0.000000 0.067256 0.029281 0.142968 331.0 0 min 10 sec
1st 19Feb2020 05Mar2020 SIR-F 4.93 0.007714 0.001322 0.136657 0.026198 0.170717 278.0 0 min 10 sec
2nd 06Mar2020 13Mar2020 SIR-F 5.06 0.022125 0.000656 0.126654 0.023830 0.093897 270.0 0 min 10 sec
3rd 14Mar2020 24Mar2020 SIR-F 2.57 0.000096 0.002956 0.063547 0.021736 0.029797 272.0 0 min 10 sec
4th 25Mar2020 02Apr2020 SIR-F 4.96 0.008054 0.000702 0.117287 0.022752 0.036893 276.0 0 min 10 sec
... ... ... ... ... ... ... ... ... ... ... ... ...
Retrospective 24th 22Nov2020 03Dec2020 SIR-F 1.29 0.000778 0.000817 0.105297 0.080746 0.014815 306.0 0 min 10 sec
25th 04Dec2020 16Dec2020 SIR-F 1.17 0.000016 0.001451 0.101517 0.085465 0.010733 377.0 0 min 10 sec
26th 17Dec2020 29Dec2020 SIR-F 1.20 0.000823 0.001372 0.100629 0.082107 0.012099 309.0 0 min 10 sec
27th 30Dec2020 31Dec2020 SIR-F 1.81 0.000399 0.001710 0.104160 0.055658 0.000481 312.0 0 min 5 sec
28th 01Jan2021 11Jun2021 SIR-F 1.00 0.000096 0.001910 0.095802 0.093858 0.360542 292.0 0 min 10 sec

72 rows × 11 columns

[60]:
# History of reproduction number
_ = snl_retro.history("Rt")
_images/usage_quick_117_0.png
[61]:
# History of Infected
_ = snl_retro.history("Infected")
_images/usage_quick_118_0.png
[62]:
# Show the representative values
snl_retro.describe()
[62]:
max(Infected) argmax(Infected) Confirmed on 11Jun2021 Infected on 11Jun2021 Fatal on 11Jun2021 2nd_Rt 4th_Rt 5th_Rt 6th_Rt 8th_Rt ... 33rd_Rt 34th_Rt 35th_Rt 36th_Rt 37th_Rt 38th_Rt 39th_Rt 40th_Rt 41st_Rt 42nd_Rt
Main 76381 14May2021 773270 34320 13253 5.06 4.96 10.91 6.38 0.17 ... 0.6 0.72 1.0 1.15 1.64 1.52 1.5 1.29 0.84 0.59
Retrospective 37443 01Jan2021 799517 35275 14598 5.17 5.08 11.26 6.06 0.18 ... - - - - - - - - - -

2 rows × 28 columns