Coronavirus propagation modeling tutorial presentation, OMLDS March 2020

Last Saturday I gave a presentation at Orlando Machine Learning and Data Science (OMLDS) Meetup.

The presentation abstract heavily borrowed descriptions and plans from a similar presentation to be given for the useR! Boston Meetup in April.

This presentation was almost entirely based on frameworks, simulations, and graphics made with Mathematica. For some parts Wolfram System Modeler was used.

Here is the (main) presentation mind-map:

(Note that mind-map’s PDF has hyperlinks.)

The presentation was given online (because of COVID-19) using Zoom. The maximum number of people registered, 100. (Many were “first timers.”) Nearly 60 showed up (and stayed throughout.)

Here is a link to the video recording: https://youtu.be/odcoi9stYuY .

References

General

[CDC1] https://www.cdc.gov/coronavirus/2019-ncov .

[WRI1] Wolfram Research, Inc. “Resources For Novel Coronavirus COVID-19”, (2020) Community.wolfram.com.

Articles, blog posts

[AA1] Anton Antonov, “Coronavirus propagation modeling considerations”, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, “Basic experiments workflow for simple epidemiological models”, (2020), SystemModeling at GitHub.

[AA3] Anton Antonov, “Scaling of Epidemiology Models with Multi-site Compartments”, (2020), SystemModeling at GitHub.

Videos

[AAv1] Anton Antonov, COVID19 Epidemic Modeling: Compartmental Models, (2020), Wolfram YouTube channel .

[AAv2] Anton Antonov, Scaling of Epidemiology Models with Multi-site Compartments, (2020), Wolfram YouTube channel .

[AAv3] Anton Antonov, Simple Economic Extension of Compartmental Epidemiological Models, (2020), Wolfram YouTube channel .

[DZv1] Diego Zviovich, Geo-spatial-temporal COVID-19 Simulations and Visualizations Over USA, (2020), Wolfram YouTube channel .

NY Times COVID-19 data visualization

Yesterday in one of the forums I frequent it was announced that New York Times has published COVID-19 data on GitHub. I decided to make a Mathematica notebook that gives data links and related code for data ingestions. (And rudimentary data analysis.)

Here is the Markdown version of the notebook: “NY Times COVID-19 data visualization”.

Here is a screenshot of the WL notebook that also links to it:

Screenshot of an interactive interface:

Histograms and Pareto principle adherence:

WirVsVirus 2020 hackathon participation

Introduction

Last weekend – 2020-03-20 ÷ 2020-03-22 – I participated in the (Germany-centric) hackathon WirVsVirus. (I friend of mine who lives in Germany asked me to team up and sign up.)

Our idea proposal was accepted, listed in the dedicated overview table (see item 806). The title of our hackathon project is:

“Geo-spatial-temporal Economic Model for COVID-19 Propagation and Management in Germany”

Nearly a dozen of people enlisted to help. (We communicated through Slack.)

13dzfagts8105

Multiple people helped with the discussion of ideas, directions where to find data, with actual data gathering, and related documented analysis. Of course, just discussing the proposed solutions was already a great help!

What was accomplished

Work plans

The following mind-map reflects pretty well what was planned and done:

15n5cjaej10q8

There is also a related org-mode file with the work plan.

Data

I obtained Germany city data with Mathematica’s build-in functions and used it to heuristically derive a traveling patterns graph, [AA1].

Here is the data:

Here is Geo-histogram of that data:

0t08vw1kjdzbc

We considered a fair amount of other data. But because of the time limitations of the hackathon we had to use only the one above.

Single-site models

During the development phase I used the model SEI2R, but since we wanted to have a “geo-spatial-temporal epidemiological economics model” I productized the implementation of SEI2HR-Econ, [AAp1].

Here are the stocks, rates, and equations of SEI2HR-Econ:

0tbp6de6zdez0

Multi-site SEI2R (SEI2HR-Econ) over a hexagonal grid graph

I managed to follow through with a large part of the work plan for the hackathon and make multi-site scaled model that “follows the money”, [AA1]. Here is a diagram that shows the travelling patterns graph and solutions at one of the nodes:

1vnygv6t7chgg

Here is (a snapshot of) an interactive interface for studying and investigating the solution results:

1pgmngb4uyuzb

For more details see the notebook [AA1]. Different parameters can be set in the “Parameters” section. Especially of interest are the quarantine related parameters: start, duration, effect on contact rates and traffic patterns.

I also put in that notebook code for exporting simulations results and programmed visualization routines in R, [AA2]. (In order other members of team to be able to explore the results.)

References

[DP1] 47_wirtschaftliche Auswirkung_Geo-spatial-temp-econ-modell, DevPost.

[WRI1] Wolfram Research, Inc., Germany city data records, (2020), SystemModeling at GitHub.

[AA1] Anton Antonov, “WirVsVirus hackathon multi-site SEI2R over a hexagonal grid graph”, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, “WirVsVirus-Hackathon in R”, (2020), SystemModeling at GitHub.

[AAp1] Anton Antonov, “Epidemiology models Mathematica package”, (2020), SystemsModeling at GitHub.

Scaling of epidemiology models with multi-site compartments

Version 1.0

Introduction

In this notebook we describe and exemplify an algorithm that allows the specification and execution geo-spatial-temporal simulations of infectious disease spreads. (Concrete implementations and examples are given.)

The assumptions of the typical compartmental epidemiological models do not apply for countries or cities that are non-uniformly populated. (For example, China, USA, Russia.) There is a need to derive epidemiological models that take into account the non-uniform distribution of populations and related traveling patterns within the area of interest.

Here is a visual aid (made with a random graph over the 30 largest cities of China):

1acjs15vamk0b
1acjs15vamk0b

In this notebook we show how to extend core, single-site epidemiological models into larger models for making spatial-temporal simulations. In the explanations and examples we use SEI2R, [AA2, AAp1], as a core epidemiological model, but other models can be adopted if they adhere to the model data structure of the package “EpidemiologyModels.m”, [AAp1].

From our experiments with we believe that the proposed multi-site extension algorithm gives a modeling ingredient that is hard emulate by other means within single-site models.

Definitions

Single-site: A geographical location (city, neighbourhood, campus) for which the assumptions of the classical compartmental epidemiological models hold.

Single site epidemiological model: A compartmental epidemiological model for a single site. Such model has a system of Ordinary Differential Equations (ODE’s) and site dependent initial conditions.

Multi-site area: An area comprised of multiple single sites with known traveling patterns between them. The area has a directed graph G with nodes that correspond to the sites and a positive matrix \text{\textit{tpm}}(G) for the traveling patterns between the sites.

Problem definition: Given (i) a single site epidemiological model M, (ii) a graph G connecting multiple sites, and (iii) a traveling patterns matrix \text{\textit{tpm}}(G) between the nodes of G derive an epidemiological model S(M,\text{\textit{tpm}}(G)) that simulates more adequately viral decease propagation over G.

Multi-Site Epidemiological Model Extension Algorithm (MSEMEA): An algorithm that derives from a given single site epidemiological model and multi-site area an epidemiological model that can be used to simulate the geo-spatial-temporal epidemics and infectious disease spreads. (The description of MSEMEA is the main message of this notebook.)

Load packages

The epidemiological models framework used in this notebook is implemented with the packages [AAp1, AAp2, AA3]; the interactive plots functions are from the package [AAp4].

Notebook structure

The section “General algorithm description” gives rationale and conceptual steps of MSEMEA.

The next two sections of the notebook follow the procedure outline using the SEI2R model as M, a simple graph with two nodes as G, and both constant and time-dependent matrices for \text{\textit{tpm}}(G).

The section “Constant traveling patterns over a grid graph” presents an important test case with a grid graph that we use to test and build confidence in MSEMEA. The sub-section “Observations” is especially of interest.

The section “Time-dependent traveling patterns over a random graph” presents a nearly “real life” application of MSEMEA using a random graph and a time dependent travelling patterns matrix.

The section “Money from lost productivity” shows how to track the money losses across the sites.

The last section “Future plans” outlines envisioned (immediate) extensions work presented in this notebook.

General algorithm description

In this section we describe a modeling approach that uses different mathematical modeling approaches for (i) the multi-site travelling patterns and (ii) the single-site disease spread interactions, and then (iii) unifies them into a common model.

Splitting and scaling

The traveling between large, densely populated cities is a very different process of the usual people mingling in those cities. The usual large, dense city mingling is assumed and used in the typical compartmental epidemiological models. It seems it is a good idea to split the two processes and derive a common model.

Assume that all journeys finish within a day. We can model the people arriving (flying in) into a city as births, and people departing a city as deaths.

Let as take a simple model like SIR or SEIR and write the equation for every site we consider. This means for every site we have the same ODE’s with site-dependent initial conditions.

Consider the traveling patterns matrix K, which is a contingency matrix derived from source-destination traveling records. (Or the adjacency matrix of a travelling patterns graph.) The matrix entry of K(i,j) tells us how many people traveled from site i to site j. We systematically change the ODE’s of the sites in following way.

Assume that site a had only travelers coming from site b and going to site b. Assume that the Total Population (TP) sizes for sites a and b are N_a and N_b respectively. Assume that only people from the Susceptible Population (SP) traveled. If the adopted single-site model is SIR, [Wk1], we take the SP equation of site a

\text{SP}_a'(t)=-\frac{\beta  \text{IP}_a(t) \text{SP}_a(t)}{N_a}-\text{SP}_a(t) \mu

and change into the equation

\text{SP}_a'(t)=-\frac{\beta  \text{IP}_a(t) \text{SP}_a(t)}{N_a}-\text{SP}_a(t) \mu -\frac{K(a,b)\text{SP}_a(t)}{N_a}+\frac{K(b,a)\text{SP}_b(t)}{N_b},

assuming that

\frac{K(a,b)\text{SP}_a(t)}{N_a}\leq N_a ,\frac{K(b,a)\text{SP}_b(t)}{N_b}\leq N_b.

Remark: In the package [AAp3] the transformations above are done with the more general and robust formula:

\min \left(\frac{K(i,j)\text{SP}_i(t)}{\text{TP}_i(t)},\text{TP}_i(t)\right).

The transformed systems of ODE’s of the sites are joined into one “big” system of ODE’s, appropriate initial conditions are set, and the “big” ODE system is solved. (The sections below show concrete examples.)

Steps of MSEMEA

MSEMEA derives a compartmental model that combines (i) a graph representation of multi-site traveling patterns with (ii) a single-site compartmental epidemiological model.

Here is a visual aid for the algorithm steps below:

0i0g3m8u08bj4
0i0g3m8u08bj4
05c2sz8hd3ryj
05c2sz8hd3ryj
  1. Get a single-site epidemiological compartmental model data structure, M.
    1. The model data structure has stocks and rates dictionaries, equations, initial conditions, and prescribed rate values; see [AA2, AAp1].
  2. Derive the site-to-site traveling patterns matrix K for the sites in the graph G.
  3. For each of node i of G make a copy of the model M and denote with M[i].
    1. In general, the models M[i], i\in G have different initial conditions.
    2. The models can also have different death rates, contact rates, etc.
  4. Combine the models M[i], i\in G into the scaled model S.
    1. Change the equations of M[i], i\in G to reflect the traveling patterns matrix K.
    2. Join the systems of ODE’s of M[i], i\in G into one system of ODE’s.
  5. Set appropriate or desired initial conditions for each of the populations in S.
  6. Solve the ODE’s of S.
  7. Visualize the results.

Precaution

Care should be taken when specifying the initial conditions of MSEMEA’s system of ODE’s (sites’ populations) and the traveling patterns matrix. For example, the simulations can “blow up” if the traveling patterns matrix values are too large. As it was indicated above, the package [AAp3] puts some safe-guards, but in our experiments with random graphs and random traveling patterns matrices occasionally we still get “wild” results.

Analogy with Large scale air-pollution modeling

There is a strong analogy between MSEMEA and Eulerian models of Large Scale Air-Pollution Modeling (LSAPM), [AA3, ZZ1].

The mathematical models of LSAPM have a “chemistry part” and an “advection-diffusion part.” It is hard to treat such mathematical model directly – different kinds of splittings are used. If we consider 2D LSAPM then we can say that we cover the modeling area with steer tank reactors, then with the chemistry component we simulate the species chemical reactions in those steer tanks, and with the advection-diffusion component we change species concentrations in the steer tanks (according to some wind patterns.)

Similarly, with MSEMEA we separated the travel of population compartments from the “standard” epidemiological modeling interaction between the population compartments.

Similarly to the so called “rotational test” used in LSAPM to evaluate numerical schemes, we derive and study the results of “grid graph tests” for MSEMEA.

Single site epidemiological model

Here is the SEI2R model from the package [AAp1]:

18y99y846b10m
18y99y846b10m

Here we endow the SEI2R model with a (prominent) ID:

0alzg909zg4h0
0alzg909zg4h0

Thus we demonstrated that we can do Step 3 of MSEMEA.

Below we use ID’s that correspond to the nodes of graphs (and are integers.)

Scaling the single-site SIR model over a small complete graph

Constant travel matrices

Assume we have two sites and the following graph and matrix describe the traveling patterns between them.

Here is the graph:

0vgm31o9drq4f
0vgm31o9drq4f

And here is the traveling patterns matrix:

0lbp0xgso2tgt
0lbp0xgso2tgt

Note that there are much more travelers from 1 to 2 than from 2 to 1.

Here we obtain the core, single-site model (as shown in the section above):

Make the multi-site compartments model with SEI2R and the two-node travel matrix using the function ToSiteCompartmentsModel of [AAp2]:

Show the unique stocks in the multi-site model:

From the symbolic form of the multi-model equations derive the specific equations with the adopted rate values:

0mjliik7acoyd
0mjliik7acoyd

Show the initial conditions:

Show the total number of equations:

Solve the system of ODE’s of the extended model:

Display the solutions for each site separately:

1o9362wmczxo6
1o9362wmczxo6

From the plots above we see that both sites start with total populations of 100000 people. Because more travelers go from 1 to 2 we see that the exposed, infected, and recovered populations are larger at 2.

Time dependent travel matrices

Instead of using constant traveling patterns matrices we can use matrices with time functions as entries. It is instructive to repeat the computations above using this matrix:

1gsoh03lixm6y
1gsoh03lixm6y

Here are the corresponding number of traveling people functions:

0qh7kbtwxyatf
0qh7kbtwxyatf

Here we scale the SIR model, solve the obtained system of ODE’s, and plot the solutions:

0trv1vnslv1rm
0trv1vnslv1rm

Note that the oscillatory nature of the temporal functions in the travelling patterns matrix are reflected in the simulation results.

Constant traveling patterns over a grid graph

In this section we do the model extension and simulation over a regular grid graph with a constant traveling patterns matrix.

Here we create a grid graph with directed edges:

0l2m5npcrlnvw
0l2m5npcrlnvw

Note that:

  • There is one directed edge between any two edge-connected nodes
  • All horizontal edges point in one direction
  • All vertical edges point in one direction
  • The edges are directed from nodes with smaller indexes to nodes with larger indexes.

Here we make a constant traveling matrix and summarize it:

0d4ocoa6gibfj
0d4ocoa6gibfj

Here we scale the SEI2R model with the grid graph constant traveling matrix:

Change the initial conditions in the following way:

  • Pick initial population size per site (same for all sites)
  • Make a constant populations vector
  • At all sites except the first one put the infected populations to zero; the first site has one severely symptomatic person
  • Set the susceptible populations to be consistent with the total and infected populations.

Solve the system of ODE’s of the scaled model:

Randomly sample the graph sites and display the solutions separately for each site in the sample:

030qbpok8qfmc
030qbpok8qfmc

Display solutions of the first and last site:

0mcxwp8l1vqlb
0mcxwp8l1vqlb

As expected from the graph structure, we can see in the first site plot that its total population is decreasing – nobody is traveling to the first site. Similarly, we can see in the last site plot that its total population is increasing – nobody leaves the last site.

Graph evolution visualizations

We can visualize the spatial-temporal evolution of model’s populations using sequences of graphs. The graphs in the sequences are copies of the multi-site graph each copy having its nodes colored according to the populations in the solutions steps.

Here is a sub-sequence for the total populations:

070ld135tkf7y
070ld135tkf7y

Here is a sub-sequence for the sum of the infected populations:

0rv8vap59g8vk
0rv8vap59g8vk

Here is a sub-sequence for the recovered population:

0dnvpy20nafpz
0dnvpy20nafpz

Here is an animation of the sum of the infected populations:

1hfd1mqh0iwk7
1hfd1mqh0iwk7

Curve shapes of the globally-aggregated solutions

Let us plot for each graph vertex the sum of the solutions of the two types of infected populations. Here is a sample of those graphs:

0jwgipeg8uznn
0jwgipeg8uznn

We can see from the plot above that at the grid vertexes we have typical SEIR curve shapes for the corresponding infected populations.

Let us evaluate the solutions for the infected populations for over all graph vertexes and sum them. Here is the corresponding “globally-aggregated” plot:

03nzzmh7yydeu
03nzzmh7yydeu

We can see that the globally aggregated plot has a very different shape than the individual vertex plots. The globally aggregated plot has more symmetric look; the individual vertex plots have much steeper gradients on their left sides.

We can conjecture that a multi-site model made by MSEMEA would capture better real life situations than any single-site model. For example, by applying MSEMEA we might be more successful in our calibration attempts for the Hubei data shown (and calibrated upon) in [AA2.].

Interactive interface

With this interactive interface we see the evolution of all populations across the graph:

0fv3dbwmah3sh
0fv3dbwmah3sh

Observations

Obviously the simulations over the described grid graph, related constant traveling patterns matrix, and constant populations have the purpose to build confidence in conceptual design of MSEMEA and its implementation.

The following observations agree with our expectations for MSEMEA’s results over the “grid graph test”.

  1. The populations plots at each site resemble the typical plots of SEI2R.
  2. The total population at the first site linearly decreases.
  3. The total population at the last site linearly increases.
  4. The plots of the total populations clearly have gradually increasing gradients from the low index value nodes to the high index value nodes.
  5. For the infected populations there is a clear wave that propagates diagonally from the low index value nodes to the high index value nodes.
    1. In the direction of the general “graph flow.“
  6. The front of the infected populations wave is much steeper (gives “higher contrast”) than the tail.
    1. This should be expected from the single-site SEI2R plots.
  7. For the recovered populations there is a clear “saturation wave” pattern that indicates that the recovered populations change from 0 to values close to the corresponding final total populations.
  8. The globally aggregated solutions might have fairly different shapes than the single-site solutions. Hence, we think that MSEMEA gives a modeling ingredient that is hard to replicate or emulate by other means in single-site models.

Time-dependent traveling patterns over a random graph

In this section we apply the model extension and simulation over a random graph with random time-dependent traveling patterns matrix.

1elzu67nqndly
1elzu67nqndly

Remark: The computations in this section work with larger random graphs; we use a small graph for more legible presentation of the workflow and results. Also, the computations should clearly demonstrate the ability to do simulations with real life data.

Derive a traveling patterns matrix with entries that are random functions:

Here is a fragment of the matrix:

17ym9q0uehfbt
17ym9q0uehfbt

Summarize and plot the matrix at t=1:

0bw9rrd64615r
0bw9rrd64615r

Here we scale the SEI2R model with the random traveling matrix:

Change the initial conditions in the following way:

  • Pick maximum population size per site
  • Derive random populations for the sites
  • At all sites except the first one put the infected populations to zero; the first site has one severely symptomatic person
  • Set the susceptible populations to be consistent with the total and infected populations.

Here solve the obtained system of ODE’s:

Here we plot the solutions:

1pzmli04rbchr
1pzmli04rbchr

Graph evolution visualizations

As in the previous section we can visualize the spatial-temporal evolution of model’s populations using sequences of graphs.

Here is a globally normalized sequence:

1v8xq8jzm6ll5
1v8xq8jzm6ll5

Here is a locally normalized (“by vertex”) sequence:

029w2jtxsyaen
029w2jtxsyaen

Money from lost productivity

The model SEI2R from [AAp1] has the stock “Money from Lost Productivity” shown as \text{MLP}(t) in the equations:

1sareh7ovkmtt
1sareh7ovkmtt

Here are MLP plots from the two-node graph model:

0cm5n3xxlewns
0cm5n3xxlewns

Here we plot the sum of the accumulated money losses:

0n6gz7j3qlq07
0n6gz7j3qlq07

Here is the corresponding “daily loss” (derivative):

0jrk2ktzeeled
0jrk2ktzeeled

Future plans

There are multiple ways to extend the presented algorithm, MSEMEA. Here are a few most immediate ones:

  1. Investigate and describe the conditions under which MSEMEA performs well, and under which it “blows up”
  2. Apply MSEMEA together with single site models that have large economics parts
  3. Do real data simulations related to the spread of COVID-19.

References

Articles, books

[Wk1] Wikipedia entry, “Compartmental models in epidemiology”.

[HH1] Herbert W. Hethcote (2000). “The Mathematics of Infectious Diseases”. SIAM Review. 42 (4): 599–653. Bibcode:2000SIAMR..42..599H. doi:10.1137/s0036144500371907.

[AA1] Anton Antonov, “Coronavirus propagation modeling considerations”, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, “Basic experiments workflow for simple epidemiological models”, (2020), SystemModeling at GitHub.

[AA3] Anton Antonov, “Air pollution modeling with gridMathematica”, (2006), Wolfram Technology Conference.

[ZZ1] Zahari Zlatev, Computer Treatment of Large Air Pollution Models. 1995. Kluwer.

Repositories, packages

[WRI1] Wolfram Research, Inc., “Epidemic Data for Novel Coronavirus COVID-19”, WolframCloud.

[AAr1] Anton Antonov, Coronavirus propagation dynamics project, (2020), SystemModeling at GitHub.

[AAp1] Anton Antonov, “Epidemiology models Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp2] Anton Antonov, “Epidemiology models modifications Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp3] Anton Antonov, “Epidemiology modeling visualization functions Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp4] Anton Antonov, “System dynamics interactive interfaces functions Mathematica package”, (2020), SystemsModeling at GitHub.

Basic experiments workflow for simple epidemiological models

Introduction

The primary purpose of this document (notebook) is to give a “stencil workflow” for simulations using the packages in the project “Coronavirus simulation dynamics”, [AAr1].

The model in this notebook – SEI2R – differs from the classical SEIR model with the following elements:

  1. Two separate infected populations: one is “severely symptomatic”, the other is “normally symptomatic”
  2. The monetary equivalent of lost productivity due to infected or died people is tracked.

Remark: We consider the coronavirus propagation models as instances of the more general System Dynamics (SD) models.

Remark: The SEI2R model is a modification of the classic epidemic model SEIR, [Wk1].

Remark: The interactive interfaces in the notebook can be used for attempts to calibrate SEI2R with real data. (For example, data for the 2019–20 coronavirus outbreak, [WRI1].)

Workflow

  1. Get one of the classical epidemiology models.
  2. Extend the equations of model if needed or desired.
  3. Set relevant initial conditions for the populations.
  4. Pick model parameters to be adjust and “play with.”
  5. Derive parametrized solutions of model’s system of equations (ODE’s or DAE’s.)
    1. Using the parameters of the previous step.
  6. Using an interactive interface experiment with different values of the parameters.
    1. In order to form “qualitative understanding.”
  7. Get real life data.
    1. Say, for the 2019-20 coronavirus outbreak.
  8. Attempt manual or automatic calibration of the model.
    1. This step will most likely require additional data transformations and programming.
    2. Only manual calibration is shown in this notebook.

Load packages of the framework

The epidemiological models framework used in this notebook is implemented in the packages [AAp1, AAp2]; the interactive plots functions are from the package [AAp3].

Getting the model code

Here we take the SEI2R model implemented in the package “EpidemiologyModels.m”, [AAp1]:

We can show a tabulated visualization of the model using the function ModelGridTableForm from [AAp1]:

0ce5juav8jq3j
0ce5juav8jq3j

Model extensions and new models

The framework implemented with the packages [AAp1, AAp2, AAp3] can be utilized using custom made data structures that follow the structure of the models in [AAp1].

Of course, we can also just extend the models from [AAp1]. In this section we show how SEI2R can be extended in two ways:

  1. By adding a birth rate added to the Susceptible Population equation (the birth rate is not included by default)
  2. By adding a new equation for the infected deceased population.

Adding births term

Here are the equations of SEI2R (from [AAp1]):

1s7f291uic6xd
1s7f291uic6xd

Here we find the position of the equation that corresponds to “Susceptible Population”:

Here we make the births term using a birth rate that is the same as the death rate:

Here we add the births term to the equations of new model

Here we display the equations of the new model:

1o2fwon3gfhel
1o2fwon3gfhel

Adding infected deceased population equation

Here we add new population, equation, and initial condition that allow for tracking the deaths because of infection:

Here is how the model looks like:

0qk5d8mdnhvu2
0qk5d8mdnhvu2

Parameters and actual simulation equations code

Here are the parameters we want to experiment with (or do calibration with):

Here we set custom rates and initial conditions:

Here is the system of ODE’s we use with to do parametrized simulations:

0dz5k6hwx6os4
0dz5k6hwx6os4

Simulation

Straightforward simulation for one year with using ParametricNDSolve :

0d6wh46looawc
0d6wh46looawc

(The advantage having parametrized solutions is that we can quickly compute simulation results with new parameter values without solving model’s system of ODE’s; see the interfaces below.)

Interactive interface

opts = {PlotRange -> All, PlotLegends -> None, PlotTheme -> "Detailed", PerformanceGoal -> "Speed", ImageSize -> 300};
lsPopulationKeys = GetPopulationSymbols[modelSI2R, __ ~~ "Population"];
lsEconKeys = {MLP};
Manipulate[
 DynamicModule[{lsPopulationPlots, lsEconPlots, lsRestPlots},
  
  lsPopulationPlots =
   ParametricSolutionsPlots[
    modelSI2R["Stocks"],
    KeyTake[aSol, lsPopulationKeys],
    {aincp, aip, spf, crisp, criap}, ndays,
    "LogPlot" -> popLogPlotQ, "Together" -> popTogetherQ, 
    "Derivatives" -> popDerivativesQ, 
    "DerivativePrefix" -> "\[CapitalDelta]", opts];
  
  lsEconPlots =
   ParametricSolutionsPlots[
    modelSI2R["Stocks"],
    KeyTake[aSol, lsEconKeys],
    {aincp, aip, spf, crisp, criap}, ndays,
    "LogPlot" -> econLogPlotQ, "Together" -> econTogetherQ, 
    "Derivatives" -> econDerivativesQ, 
    "DerivativePrefix" -> "\[CapitalDelta]", opts];
  
  lsRestPlots =
   ParametricSolutionsPlots[
    modelSI2R["Stocks"],
    KeyDrop[aSol, Join[lsPopulationKeys, lsEconKeys]],
    {aincp, aip, spf, crisp, criap}, ndays,
    "LogPlot" -> econLogPlotQ, "Together" -> econTogetherQ, 
    "Derivatives" -> econDerivativesQ, 
    "DerivativePrefix" -> "\[CapitalDelta]", opts];
  
  Multicolumn[Join[lsPopulationPlots, lsEconPlots, lsRestPlots], 
   nPlotColumns, Dividers -> All, 
   FrameStyle -> GrayLevel[0.8]]
  ],
 {{aincp, 12., "Average incubation period (days)"}, 1, 60., 1, Appearance -> {"Open"}},
 {{aip, 21., "Average infectious period (days)"}, 1, 100., 1, Appearance -> {"Open"}},
 {{spf, 0.2, "Severely symptomatic population fraction"}, 0, 1, 0.025, Appearance -> {"Open"}},
 {{crisp, 6, "Contact rate of the infected severely symptomatic population"}, 0, 30, 0.1, Appearance -> {"Open"}},
 {{criap, 3, "Contact rate of the infected normally symptomatic population"}, 0, 30, 0.1, Appearance -> {"Open"}},
 {{ndays, 90, "Number of days"}, 1, 365, 1, Appearance -> {"Open"}},
 {{popTogetherQ, True, "Plot populations together"}, {False, True}},
 {{popDerivativesQ, False, "Plot populations derivatives"}, {False, True}},
 {{popLogPlotQ, False, "LogPlot populations"}, {False, True}},
 {{econTogetherQ, False, "Plot economics functions together"}, {False, True}},
 {{econDerivativesQ, False, "Plot economics functions derivatives"}, {False, True}},
 {{econLogPlotQ, False, "LogPlot economics functions"}, {False, True}},
 {{nPlotColumns, 1, "Number of plot columns"}, Range[5]},
 ControlPlacement -> Left, ContinuousAction -> False]
0uhcbh5jg8g3a
0uhcbh5jg8g3a

Calibration over real data

It is important to calibrate these kind of models with real data, or at least to give a serious attempt to such a calibration. If the calibration is “too hard” or “impossible” that would indicate that the model is not that adequate. (If adequate at all.)

The calibration efforts can be (semi-)automated using special model-to-data goodness of fit functions and a minimization procedure. (See for example, [AA2].)

In this section we just attempt to calibrate SEI2R over real data taken from [WRI1] using a specialized interactive interface.

Real data

Here is COVID-19 data taken from [WRI1] for the Chinese province Hubei:

The total population in Hubei is

1kt1ikvs8tzqt
1kt1ikvs8tzqt
1cpkt5fvgh8hu
1cpkt5fvgh8hu

But we have to use a fraction of that population in order to produce good fits. That can be justified with the conjecture that the citizens of Hubei are spread out and it is mostly one city (Wuhan) where the outbreak is.

The real data have to be padded with a certain number of 0’s related to the infectious and incubation periods in order to make good fits. Such padding is easy to justify: if we observe recovered people that means that they passed through the incubation and infectious periods.

Calibration interactive interface

In this interface we put the Infected Severely Symptomatic Population (ISSP) to zero. That way it is easier to compare the real data with the simulated results (and pick parameter values that give close fits.) Also note that since SEI2R is simple in this interface the system is always solved.

opts = {PlotRange -> All, PlotLegends -> None, PlotTheme -> "Detailed", PerformanceGoal -> "Speed", ImageSize -> 300};
Manipulate[
 DynamicModule[{modelSI2R = modelSI2R, lsActualEquations, aSol, 
   lsPopulationPlots, lsEconPlots, lsRestPlots},
  
  modelSI2R = SetRateRules[modelSI2R, <|TP[t] -> population|>];
  modelSI2R = 
   SetInitialConditions[
    modelSI2R, <|SP[0] -> population - 1, ISSP[0] -> 0, 
     INSP[0] -> 1|>];
  lsActualEquations = 
   Join[modelSI2R["Equations"] //. 
     KeyDrop[modelSI2R["RateRules"], lsFocusParams], 
    modelSI2R["InitialConditions"]];
  aSol =
   Association@Flatten@
     ParametricNDSolve[
      lsActualEquations, {SP, EP, INSP, RP, IDP}, {t, 0, 365}, 
      lsFocusParams];
  
  lsPopulationPlots =
   ParametricSolutionsPlots[
    modelSI2R["Stocks"],
    KeyTake[aSol, GetPopulationSymbols[modelSI2R, __ ~~ "Population"]],
    {aincp, aip, 0, criap, criap}, ndays, "Together" -> True, 
    opts];
  
  Show[lsPopulationPlots[[1]], 
   ListPlot[
    PadRealData[aRealData, Round[aincp + padOffset], 
     Round[aip + padOffset]], PlotStyle -> {Blue, Black, Red}]]
  ],
 {{population, 58160000/600, "Population"}, 58160000/1000, 58160000, 10000, Appearance -> {"Open"}},
 {{padOffset, 0, "real data padding offset"}, -100, 100, 1, Appearance -> {"Open"}},
 {{aincp, 6, "Average incubation period (days)"}, 1, 60, 1, Appearance -> {"Open"}},
 {{aip, 32, "Average infectious period (days)"}, 1, 100, 1, Appearance -> {"Open"}},
 {{criap, 0.8, "Contact rate of the infected normally symptomatic population"}, 0, 30, 0.1, Appearance -> {"Open"}},
 {{ndays, 90, "Number of days"}, 1, 365, 1, Appearance -> {"Open"}},
 ControlPlacement -> Left, ContinuousAction -> False]
0s4dnliwjni2v
0s4dnliwjni2v

Maybe good enough parameters

1v43idv1zv24j
1v43idv1zv24j

Basic reproduction number:

(*25.5966*)
0upbzla7bc2ok
0upbzla7bc2ok

Basic reproduction number:

(*59.7934*)

References

Articles

[Wk1] Wikipedia entry, “Compartmental models in epidemiology”.

[HH1] Herbert W. Hethcote (2000). “The Mathematics of Infectious Diseases”. SIAM Review. 42 (4): 599–653. Bibcode:2000SIAMR..42..599H. doi:10.1137/s0036144500371907.

[AA1] Anton Antonov, “Coronavirus propagation modeling considerations”, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, Answer of “Model calibration with phase space data”, (2019), Mathematica StackExchage.

Repositories & packages

[WRI1] Wolfram Research, Inc., “Epidemic Data for Novel Coronavirus COVID-19”, WolframCloud.

[AAr1] Anton Antonov, Coronavirus propagation dynamics project, (2020), SystemModeling at GitHub.

[AAp1] Anton Antonov, “Epidemiology models Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp2] Anton Antonov, “Epidemiology models modifications Mathematica package”, (2020), SystemsModeling at GitHub.

[AAp3] Anton Antonov, “System dynamics interactive interfaces functions Mathematica package”, (2020), SystemsModeling at GitHub.

Conference abstracts similarities

Introduction

In this MathematicaVsR project we discuss and exemplify finding and analyzing similarities between texts using Latent Semantic Analysis (LSA). Both Mathematica and R codes are provided.

The LSA workflows are constructed and executed with the software monads LSAMon-WL, [AA1, AAp1], and LSAMon-R, [AAp2].

The illustrating examples are based on conference abstracts from rstudio::conf and Wolfram Technology Conference (WTC), [AAd1, AAd2]. Since the number of rstudio::conf abstracts is small and since rstudio::conf 2020 is about to start at the time of preparing this project we focus on words and texts from RStudio’s ecosystem of packages and presentations.

Statistical thesaurus for words from RStudio’s ecosystem

Consider the focus words:

{"cloud","rstudio","package","tidyverse","dplyr","analyze","python","ggplot2","markdown","sql"}

Here is a statistical thesaurus for those words:

0az70qt8noeqf
0az70qt8noeqf

Remark: Note that the computed thesaurus entries seem fairly “R-flavored.”

Similarity analysis diagrams

As expected the abstracts from rstudio::conf tend to cluster closely – note the square formed top-left in the plot of a similarity matrix based on extracted topics:

1d5a83m8cghew
1d5a83m8cghew

Here is a similarity graph based on the matrix above:

09y26s6kr3bv9
09y26s6kr3bv9

Here is a clustering (by “graph communities”) of the sub-graph highlighted in the plot above:

0rba3xgoknkwi
0rba3xgoknkwi

Notebooks

Comparison observations

LSA pipelines specifications

The packages LSAMon-WL, [AAp1], and LSAMon-R, [AAp2], make the comparison easy – the codes of the specified workflows are nearly identical.

Here is the Mathematica code:

lsaObj =
  LSAMonUnit[aDesriptions]⟹
   LSAMonMakeDocumentTermMatrix[{}, Automatic]⟹
   LSAMonEchoDocumentTermMatrixStatistics⟹
   LSAMonApplyTermWeightFunctions["IDF", "TermFrequency", "Cosine"]⟹
   LSAMonExtractTopics["NumberOfTopics" -> 36, "MinNumberOfDocumentsPerTerm" -> 2, Method -> "ICA", MaxSteps -> 200]⟹
   LSAMonEchoTopicsTable["NumberOfTableColumns" -> 6];

Here is the R code:

lsaObj <- 
  LSAMonUnit(lsDescriptions) %>% 
  LSAMonMakeDocumentTermMatrix( stemWordsQ = FALSE, stopWords = stopwords::stopwords() ) %>% 
  LSAMonApplyTermWeightFunctions( "IDF", "TermFrequency", "Cosine" ) 
  LSAMonExtractTopics( numberOfTopics = 36, minNumberOfDocumentsPerTerm = 5, method = "NNMF", maxSteps = 20, profilingQ = FALSE ) %>% 
  LSAMonEchoTopicsTable( numberOfTableColumns = 6, wideFormQ = TRUE ) 

Graphs and graphics

Mathematica’s built-in graph functions make the exploration of the similarities much easier. (Than using R.)

Mathematica’s matrix plots provide more control and are more readily informative.

Sparse matrix objects with named rows and columns

R’s built-in sparse matrices with named rows and columns are great. LSAMon-WL utilizes a similar, specially implemented sparse matrix object, see [AA1, AAp3].

References

Articles

[AA1] Anton Antonov, A monad for Latent Semantic Analysis workflows, (2019), MathematicaForPrediction at GitHub.

[AA2] Anton Antonov, Text similarities through bags of words, (2020), SimplifiedMachineLearningWorkflows-book at GitHub.

Data

[AAd1] Anton Antonov, RStudio::conf-2019-abstracts.csv, (2020), SimplifiedMachineLearningWorkflows-book at GitHub.

[AAd2] Anton Antonov, Wolfram-Technology-Conference-2016-to-2019-abstracts.csv, (2020), SimplifiedMachineLearningWorkflows-book at GitHub.

Packages

[AAp1] Anton Antonov, Monadic Latent Semantic Analysis Mathematica package, (2017), MathematicaForPrediction at GitHub.

[AAp2] Anton Antonov, Latent Semantic Analysis Monad R package, (2019), R-packages at GitHub.

[AAp3] Anton Antonov, SSparseMatrix Mathematica package, (2018), MathematicaForPrediction at GitHub.

Pets licensing data analysis

Introduction

This notebook / document provides ground data analysis used to make or confirm certain modeling conjectures and assumptions of a Pets Retail Dynamics Model (PRDM), [AA1]. Seattle pets licensing data is used, [SOD2].

We want to provide answers to the following questions.

  • Does the Pareto principle manifests for pets breeds?

  • Does the Pareto principle manifests for ZIP codes?

  • Is there an upward trend for becoming a pet owner?

All three questions have positive answers, assuming the retrieved data, [SOD2], is representative. See the last section for an additional discussion.

We also discuss pet adoption simulations that are done using Quantile Regression, [AA2, AAp1].

This notebook/document is part of the SystemsModeling at GitHub project “Pets retail dynamics”, [AA1].

Data

The pet licensing data was taken from this page: “Seattle Pet Licenses”, https://data.seattle.gov/Community/Seattle-Pet-Licenses/jguv-t9rb/data.

The ZIP code coordinates data was taken from a GitHub repository,
“US Zip Codes from 2013 Government Data”, https://gist.github.com/erichurst/7882666.

Animal licenses

image-3281001a-2f3d-4a8e-87b9-dc8a8b9803b3
image-3281001a-2f3d-4a8e-87b9-dc8a8b9803b3

Convert “Licence Issue Date” values into DateObjects.

Summary

image-49aecba4-2b43-40d7-87ba-15ceb848898d
image-49aecba4-2b43-40d7-87ba-15ceb848898d

Keep dogs and cats only

Since the number of animals that are not cats or dogs is very small we remove them from the data in order to produce more concise statistics.

ZIP code geo-coordinates

Summary

image-572ef441-b14e-438d-b5b7-85f244aa1857
image-572ef441-b14e-438d-b5b7-85f244aa1857
image-c0d4f154-ee22-457f-8a36-715b77c92e08
image-c0d4f154-ee22-457f-8a36-715b77c92e08

Pareto principle adherence

In this section we apply the Pareto principle statistic in order to see does the Pareto principle manifests over the different columns of the pet licensing data.

Breeds

We see a typical Pareto principle adherence for both dog breeds and cat breeds. For example, 20% of the dog breeds correspond to 80% of all registered dogs.

Note that the number of unique cat breeds is 4 times smaller than the number of unique dog breeds.

image-d1bac8f8-fe6c-42c0-8d52-45ed21ab6cc2
image-d1bac8f8-fe6c-42c0-8d52-45ed21ab6cc2
image-3c320985-1ed4-4d11-b983-29f87d4cdc7c
image-3c320985-1ed4-4d11-b983-29f87d4cdc7c

Animal names

We see a typical Pareto principle adherence for the frequencies of the pet names. For dogs, 10% of the unique names correspond to ~65% of the pets.

image-cb6368b6-b735-4f77-a3dd-bcb0be60f28e
image-cb6368b6-b735-4f77-a3dd-bcb0be60f28e
image-bbcac6bb-5247-400c-a093-f3002206b5cf
image-bbcac6bb-5247-400c-a093-f3002206b5cf

Zip codes

We see typical – even exaggerated – manifestation of the Pareto principle over ZIP codes of the registered pets.

image-72cae8dd-d342-4c90-a11d-11607545133e
image-72cae8dd-d342-4c90-a11d-11607545133e

Geo-distribution

In this section we visualize the pets licensing geo-distribution with geo-histograms.

Both cats and dogs

image-94ae1316-ada2-4195-b2fc-6864ff1fd642
image-94ae1316-ada2-4195-b2fc-6864ff1fd642

Dogs and cats separately

image-836dff19-7000-45e0-b0a4-1f3fe4a066c9
image-836dff19-7000-45e0-b0a4-1f3fe4a066c9

Pet stores

In this subsection we show the distribution of pet stores (in Seattle).

It is better instead of image retrieval to show corresponding geo-markers in the geo-histograms above. (This is not considered that important in the first version of this notebook/document.)

image-836dff19-7000-45e0-b0a4-1f3fe4a066c9
image-836dff19-7000-45e0-b0a4-1f3fe4a066c9

Time series

In this section we visualize the time series corresponding to the pet registrations.

Time series objects

Here we make time series objects:

image-49ae54cb-0644-427e-a015-0392284aaaa7
image-49ae54cb-0644-427e-a015-0392284aaaa7

Time series plots of all registrations

Here are time series plots corresponding to all registrations:

image-02632be6-ab52-41b8-959a-e200641fdd8f
image-02632be6-ab52-41b8-959a-e200641fdd8f

Time series plots of most recent registrations

It is an interesting question why the number of registrations is much higher in volume and frequency in the years 2018 and later.

image-85ebeab1-cad5-4fe3-bd5d-c7c8c94a753e
image-85ebeab1-cad5-4fe3-bd5d-c7c8c94a753e

Upward trend

Here we apply both Linear Regression and Quantile Regression:

image-6df4d9d2-e48a-4d63-885c-6ed5112c0f15
image-6df4d9d2-e48a-4d63-885c-6ed5112c0f15

We can see that there is clear upward trend for both dogs and cats.

Quantile regression application

In this section we investigate the possibility to simulate the pet adoption rate. We plan to use simulations of the pet adoption rate in PRDM.

We do that using the software monad QRMon, [AAp1]. A list of steps follows.

  • Split the time series into windows corresponding to the years 2018 and 2019.

  • Find the difference between the two years.

  • Apply Quantile Regression to the difference using a reasonable grid of probabilities.

  • Simulate the difference.

  • Add the simulated difference to year 2019.

Simulation

In this sub-section we simulate the differences between the time series for 2018 and 2019, then we add the simulated difference to the time series of the year 2019.

image-8f9e3af0-46b7-4417-bd1e-3201c1134f34
image-8f9e3af0-46b7-4417-bd1e-3201c1134f34
image-30b836dc-f166-4f21-9c0b-9cca922058e6
image-30b836dc-f166-4f21-9c0b-9cca922058e6
image-65e4d1bf-dfff-4073-88a0-63177eeed1b5
image-65e4d1bf-dfff-4073-88a0-63177eeed1b5
image-6d107cad-6fef-46c8-92a8-59ea78b5039f
image-6d107cad-6fef-46c8-92a8-59ea78b5039f
image-d0d517e0-925b-486c-88fd-287cfe02e799
image-d0d517e0-925b-486c-88fd-287cfe02e799

Take the simulated time series difference:

Add the simulated time series difference to year 2019, clip the values less than zero, shift the result to 2020:

image-2a29feca-73b8-4fce-8051-145d74ec499c
image-2a29feca-73b8-4fce-8051-145d74ec499c

Plot all years together

image-793f146a-07f9-455f-9bc7-2ef7d7897691
image-793f146a-07f9-455f-9bc7-2ef7d7897691

Discussion

This section has subsections that correspond to additional discussion questions. Not all questions are answered, the plan is to progressively answer the questions with the subsequent versions of the this notebook / document.

□ Too few pets

The number of registered pets seems too few. Seattle is a large city with more than 600000 citizens; approximately 50% of the USA households have dogs; hence the registered pets are too few (~50000).

□ Why too few pets?

Seattle is a high tech city and its citizens are too busy to have pets?

Most people do not register their pets? (Very unlikely if they have used veterinary services.)

Incomplete data?

□ Registration rates

Why the number of registrations is much higher in volume and frequency in the years 2018 and later?

□ Adoption rates

Can we tell apart the adoption rates of pet-less people and people who already have pets?

Preliminary definitions

References

[AA1] Anton Antonov, Pets retail dynamics project, (2020), SystemModeling at GitHub.

[AA2] Anton Antonov, A monad for Quantile Regression workflows, (2018), MathematicaForPrediction at WordPress.

[AAp1] Anton Antonov, Monadic Quantile Regression Mathematica package, (2018), MathematicaForPrediction at GitHub.

[SOD1] Seattle Open Data, “Seattle Pet Licenses”, https://data.seattle.gov/Community/Seattle-Pet-Licenses/jguv-t9rb/data .