Message Boards

GROUPS:

Staff Picks Biological Sciences Data Science Mathematics Medical Sciences Graphics and Visualization Wolfram Language Modeling Statistics and Probability Events & Media

Robert Rimmer

[Notebook] Updated: coronavirus logistic growth model: China

Robert Rimmer

Posted 7 months ago

3125 Views

7 Replies

11 Total Likes

Follow this post

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus

enter image description here

On 29 February 2020, the China National Health Commission reported 423 new cases of coronavirus for Hubei and only 4 cases for the rest of the country. The logistic function fit to the data from the Johns Hopkins CSSE, have been showing for a while that the epidemic should end outside of the Hubei province in the first week of March. The recent report seems to confirm this prediction.

This notebook contains some updated and new functions that were not in my previous posts.

Data Source

https://gisanddata.maps.arcgis.com/apps/opsdashboard/index.html#/bda7594740fd40299423467b48e9ecf6

The latest data gets updated on dashboard above before it is loaded to GitHub, which is the source for the notebook.

https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_ 19_data/csse_covid_ 19_time _series

Procedures

Clear[modelFitDisplay];modelFitDisplay[NLM_,data_,opts___]:=Module[{bands,p0,p1},bands[t_]=NLM["SinglePredictionBands",ConfidenceLevel0.95];p0=Plot[{NLM[t],bands[t]},{t,1,data[[-1]][[1]]+5},opts,Filling{2{1}},PlotStyle{Darker@Red,Darker@Orange,Darker@Orange},GridLinesAutomatic,PlotRangeAll,opts];p1=ListPlot[Tooltip[#,#]&/@data,PlotRangeAll];Show[p0,p1]];

In[]:=

Clear[derivativeDisplay];derivativeDisplay[NLM_,data_,opts___]:=Module[{daynum,end,peak},daynum=data[[-1]][[1]];peak=t0/.NLM["BestFitParameters"];end=x/.FindRoot[NLM'[x]1,{x,peak,peak+200},(*forcetorighttail*)Method"Brent"][[1]];Print["Last case: ",DateString[DatePlus[{2020,1,21},Ceiling@end],"DateShort"]];Show[Plot[{NLM'[x],NLM''[x]},{x,data[[1]][[1]]-4,daynum+5},opts,GridLinesAutomatic,PlotRangeAll,PlotLegends{"First Derivative","Second Derivative"}],ListPlot[Tooltip[#,#]&/@{{daynum,NLM'[daynum]}},PlotStyle{Red,PointSize[0.02]}]]];

In[]:=

Clear[dataD1];dataD1[data_,{k_,L_}]:=Table[{data[[i]][[1]],kdata[[i]][[2]](1-data[[i]][[2]]/L)},{i,Length[data]}];

In[]:=

logistic=L/(1+Exp[-k(t-t0)]);

In[]:=

Clear[logisticDEFit];logisticDEFit[data_,graphic_:True,optsInterpolation___]:=Module{dataFn,dataDE,dataDEPlot,deModel,deNLM,bands},dataFn=Interpolation[data,optsInterpolation,Method"Spline"];dataDE={dataFn[#],dataFn'[#]}&/@data[[All,1]];dataDEPlot=ListPlot[dataDE,PlotStyleDarker@Red,GridLinesAutomatic,PlotRangeAll];deModel=kx1-

;deNLM=NonlinearModelFit[dataDE,deModel,{k,L},x];If[graphic,bands[L_]=deNLM["SinglePredictionBands",ConfidenceLevel0.95];Print@Show[Plot[{deNLM[x],bands[x]},{x,0,L/.deNLM@"BestFitParameters"},GridLinesAutomatic,Filling{2->{1}},FrameTrue,FrameLabel{"Cases","New Cases Per Day"},PlotRangeAll,PlotLabel"Logistic Differential Equation"],dataDEPlot]];Flatten@{deNLM@"BestFitParameters",FindRoot[dataFn[t0]L/2/.deNLM["BestFitParameters"],{t0,data[[-1]][[1]]}]};

In[]:=

Import Data

deaths=Import["https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Deaths.csv"];deaths[[1]][[-1]]

In[]:=

2/28/20

Out[]=

deathData=AssociationThread[First@deaths,#]&/@Rest[deaths]//Dataset;

In[]:=

confirmed=Import["https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv"];confirmed[[1]][[-1]]

In[]:=

2/28/20

Out[]=

Setup Datasets

caseData=AssociationThread[First@confirmed,#]&/@Rest[confirmed]//Dataset;

In[]:=

chinaCaseData=caseData@Select[#"Country/Region""Mainland China"&];chinaDeathData=deathData@Select[#"Country/Region""Mainland China"&];

In[]:=

keys=Keys[First@Normal@caseData];Keys[First@Normal@deathData];

In[]:=

deathCases=Transpose[{Range[Length@keys-4],Values@Normal[chinaDeathData[Total,Drop[keys,4]]]}];confirmedCases=Transpose[{Range[Length@keys-4],Values@Normal[chinaCaseData[Total,Drop[keys,4]]]}];

In[]:=

notHubeiData=chinaCaseData@Select[#"Province/State"≠"Hubei"&];notHubeiCases=Transpose[{Range[Length@keys-4],Values@Normal[notHubeiData[Total,Drop[keys,4]]]}];

In[]:=

hubeiData=chinaCaseData@Select[#"Province/State""Hubei"&];hubeiCases=Transpose[{Range[Length@keys-4],Values@Normal[hubeiData[All,Drop[keys,4]]][[1]]}];

In[]:=

China Death Data

deathsNLM=NonlinearModelFit[deathCases,logistic,{{k,0.3},{L,3000},{t0,25}},t];deathsNLM@"BestFitParameters"

In[]:=

{k0.172041,L3044.88,t024.2041}

Out[]=

deathsNLM@"ParameterConfidenceIntervalTable"

In[]:=

	Estimate	Standard Error	Confidence Interval
k	0.172041	0.00394588	{0.16403,0.180052}
L	3044.88	41.6121	{2960.4,3129.36}
t0	24.2041	0.226153	{23.745,24.6632}

Out[]=

The graphic below shows the first and second derivative plots of the logistic function. The red dot is at the current day number on the first derivative curve. The green dots are the points calculated by the logistic differential equation, using each data point. The prediction for the last case date is printed (probably not accurate until there are at least three weeks of new data).

Last case: Sun 22 Mar 2020

	First Derivative
	Second Derivative

Out[]=

China Clinical Cases

casesNLM=NonlinearModelFit[confirmedCases,logistic,{{k,0.3},{L,80000},{t0,25}},t];casesNLM@"BestFitParameters"

In[]:=

{k0.221968,L80902.5,t018.7671}

Out[]=

casesNLM@"ParameterConfidenceIntervalTable"

In[]:=

	Estimate	Standard Error	Confidence Interval
k	0.221968	0.0103505	{0.200955,0.24298}
L	80902.5	1230.14	{78405.2,83399.8}
t0	18.7671	0.265818	{18.2274,19.3067}

Out[]=

Last case: Tue 24 Mar 2020

	First Derivative
	Second Derivative

Out[]=

China Cases Outside Hubei

notHubeiNLM=NonlinearModelFit[notHubeiCases,logistic,{{k,0.3},{L,15000},{t0,25}},t];notHubeiNLM@"BestFitParameters"

In[]:=

{k0.25761,L12801.6,t013.4586}

Out[]=

notHubeiNLM@"ParameterConfidenceIntervalTable"

In[]:=

	Estimate	Standard Error	Confidence Interval
k	0.25761	0.00547199	{0.246502,0.268719}
L	12801.6	62.0486	{12675.7,12927.6}
t0	13.4586	0.0955942	{13.2645,13.6526}

Out[]=

Last case: Fri 6 Mar 2020

	First Derivative
	Second Derivative

Out[]=

China Hubei Cases

hubeiNLM=NonlinearModelFit[hubeiCases,logistic,{{k,0.3},{L,15000},{t0,25}},t];hubeiNLM@"BestFitParameters"

In[]:=

{k0.235356,L67672.4,t019.6809}

Out[]=

hubeiNLM@"ParameterConfidenceIntervalTable"

In[]:=

	Estimate	Standard Error	Confidence Interval
k	0.235356	0.0127775	{0.209417,0.261296}
L	67672.4	1189.32	{65257.9,70086.8}
t0	19.6809	0.292337	{19.0874,20.2744}

Out[]=

Last case: Sun 22 Mar 2020

	First Derivative
	Second Derivative

Out[]=

Fitting to the Logistic Differential Equation

These graphics show a method of fitting to the equation for the logistic first derivative.

f'[t]=kx1-

The value of x is the number of cumulative cases. This is a quadratic equation with roots at x = 0 and x = L, and peak at x = L/2. The first derivative is obtained from the data by making an interpolation function (Method  “Spline” is default) from the data then using the function to calculate first derivative at each point. The 95% confidence intervals are shown for the quadratic fit. Clearly the data sets show quite a lot of irregularity by this method. The smoother data set, notHubei, shows the later points becoming more accurate. The parameters are similar to the fit to the logistic function.

logisticDEFit[deathCases]

In[]:=

{k0.167235,L3097.5,t024.2121}

Out[]=

logisticDEFit[confirmedCases]

In[]:=

{k0.221009,L80819.9,t019.2458}

Out[]=

logisticDEFit[notHubeiCases]

In[]:=

{k0.252988,L12850.2,t045.7223}

Out[]=

logisticDEFit[hubeiCases]

In[]:=

{k0.237076,L67516.9,t022.0443}

Out[]=

POSTED BY: Robert Rimmer

Answer

7 Replies

Sort By:

Mads Bahrami

Mads Bahrami, Wolfram Research Inc

Posted 7 months ago

Thanks a lot Robert! I really enjoy reading it. I wonder why you have used interpolation to get dataDE in your function logisticDEFit? For example, look at this notebook where I compare my suggestion with your interpolation data:

POSTED BY: Mads Bahrami

Answer

Robert Rimmer

Posted 7 months ago

POSTED BY: Robert Rimmer

Answer

Mads Bahrami

Mads Bahrami, Wolfram Research Inc

Posted 7 months ago

Thanks for your input! Well-explained!

POSTED BY: Mads Bahrami

Answer

Moderation Team

Moderation Team, WOLFRAM

Posted 7 months ago

-- you have earned *Featured Contributor Badge* Your exceptional post has been selected for our editorial column *Staff Picks* http://wolfr.am/StaffPicks and Your Profile is now distinguished by a *Featured Contributor Badge* and is displayed on the Featured Contributor Board. Thank you, keep it coming, and consider contributing to the The Notebook Archive!

POSTED BY: Moderation Team

Answer

Robert Rimmer

Posted 7 months ago

I posted a notebook tying the SEIR model to the Logistic model: LogisticFromSEIR

POSTED BY: Robert Rimmer

Answer

Updating Name

Posted 7 months ago

Hi, I really appreciate the job you did. I use to use statistics in my professional activity, and I'm quite familiar with problems concerning the "bayesian calibration of computer code outputs". Since a few weeks I'm analysing the data about the CoVid19 outbreak and one of my main question is "how can we infer the value of the mortality rate" ? This is not a trivial problem and I think that bayesian inference (I don't know if you are familiar with this topic?) could enable estimating the distribution of "credible" mortality rates. Do you have some knowledge about works concerning this approach? df

POSTED BY: Updating Name

Answer

Robert Rimmer

Posted 7 months ago

Right now I don't think the problem can be solved. Deaths are easy to count and because they were ill first most all were tested, but the number of infected people is determined by testing, which has varied widely by location and indication. Until there has been surveillance testing in a large community to know the incidence and prevalence of infection, I don't think you can estimate probabilities needed for Bayesian inference.

POSTED BY: Robert Rimmer

Answer

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback