File Exchange

image thumbnail

fitVirusCV19v3 (COVID-19 SIR Model)

version 3.4.8 (3.57 MB) by Joshua McGee
Estimation of coronavirus COVID-19 epidemic evaluation by the SIR model, Code receives and plots most recent data from HDX.

475 Downloads

Updated 03 May 2020

View License

This code will retrieve recent data from the HDX and attempt to fit it to a SIR model using structure created by Milan Batista (https://www.mathworks.com/matlabcentral/fileexchange/74658-fitviruscovid19?s_tid=prof_contriblnk). The model is data-driven, so its forecast is as good as data are. Also, it is assumed that the model is a reasonable description of the one-stage epidemic.

Results are saved in structure res (see function fiVirusCV19 header). Optionally the results may be printed by:

fitVirusCV19v3(“Italy",'prn','on')

The plot may be enabled or disabled via (default is on):
fitVirusCV19v3(“Italy",'plt','off')

To plot the growth rate on the figure use (def value is 2)
fitVirusCV19v3("Italy",'nsp',3)

To plot the R0 as a function of time
fitVirusCV19v3("Italy",'r0analysis','on')

A more detailed description can be found in:
https://www.researchgate.net/publication/339311383_Estimation_of_the_final_size_of_the_coronavirus_epidemic_by_the_SIR_model
Examples can be found in:
https://www.researchgate.net/publication/339912313_Forecasting_of_final_COVID-19_epidemic_size_200318

Cite As

Joshua McGee (2020). fitVirusCV19v3 (COVID-19 SIR Model) (https://www.mathworks.com/matlabcentral/fileexchange/74676-fitviruscv19v3-covid-19-sir-model), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (98)

Autumn Lee

Autumn Lee

These are the original output from Milan Batista SIR model for Singapore latest data , I don't think your current programme is able to generate result for Singapore , please revisit all your coding and assumption again , I have tested a few and noticed your final result appeared to differ from the original Milan Bastista computation for most of the countries I have tested.

fitVirusCV19(@getDataSingapore,'plt','on','prn','on','data','on');

Epidemic modeling by susceptible-infected-recovered (SIR) model
Country Singapore
Day 118
Estimated the SIR model parameters
Contact frequency (beta) 0.244 (1/day)
Removal frequency (gamma) 0.123 (1/day)
Population size (N) 42561
Initial number of cases (I0) 0
Basic rep. number (R0) 1.989
Reproduction number (R ) 0.59
End reproduction number 0.41
Time between contacts (Tc) 4.1 (day)
Infectious period (Tr) 8.1 (day)
Final state
Final number of cases 33779
Final number of susceptibles 8781
Daily forecast for 25-May-2020
Total NaN
Increase NaN
Estimated logistic model parameters
Epidemic size (K) 28265 (cases)
Epidemic rate (r) 0.121429 (1/day)
Initial doubling time 5.7 (day)
Estimated duration (days)
Turning day 93
Acceleration phase 17 (days)
Deceleration phase 20 (days)
Total growth phase duration 37 (days)
Total epidemic duration 230 (days)
Estimated datums
Outbreak 28-Jan-2020
Start of acceleration 13-Apr-2020
Turning point 30-Apr-2020
Start of steady growth 21-May-2020
Start of ending phase 09-Jun-2020
End of epidemic (5 cases) 23-Aug-2020
End of epidemic (1 case) 14-Sep-2020
Statistics (total cases)
Number of observations 118
Degrees of freedom 114
Root Mean Squared Error 484.562
R-Squared 0.998
Adjusted R-Squared 0.998
F-statistics vs. zero model 15782.6
p-value 4.2134e-149
Statistics (daily new cases)
Number of observations 117
Degrees of freedom 113
Root Mean Squared Error 180.421
R-Squared 0.758
Adjusted R-Squared 0.75
F-statistics vs. zero model 119.252
p-value 7.22686e-35
Method
Total cases weight 0.5
Infection rate weight 0.5
Objective function value 3545.37
Exit condition (1=OK)

Autumn Lee

@Yafia Radouane I think the reason why R_0 never decrease is may be because the model compute R using the running total cumulative 4 parameters logistics Curve approximate models without discounting the total actual number of recover person, I am not sure if this is in fact default SIR model of Milan Batista SIR =S+I+R where R is recovered,RN = beta/gamma -- basic reproduction number (ratio) , only Joshua can clarity as there is the no cross validation against the actual recovery data reported by each country. In my empirical studies , 5 parameters logistics Curve seem to be a better fit to majority of the actual SARS-CoV-2 cases such as L5P from https://www.mathworks.com/matlabcentral/fileexchange/38043.

Hi Joshua and everyone on this blog;
I think this is a problem in computing R_0, because even if the number of infected cases decreases the value of R_0 still bigger than 1 which not normal. Can some one can explain this situation?
Dear Joshua, please can give me an answer???
Many thanks

I am getting the following error:

Error using fitVirusCV19v3 (line 174)
Table variable names must be character vectors.

Please can anybody help!!

Thank you for this work. it works the code well. It's awesome.

Autumn Lee

Just trying to clarify on gamma which is supposed to be recovery rate, I can not find any actual recovery data being applied to cross validate gamma computed? Could you be able to share info on this especially gamma in your computation is termed removal rate not recovery rate?

Autumn Lee

Testing on USA , Taiwan, Singapore y yield errors, could this due to rigid 4 parameters logistics curve bring use? Could a 5 parameters Richard logistics curve be added for better model fits?

Autumn Lee

Just for sharing , the result.txt snd deathresult.txt. Could also be obtain from https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data/csse_covid_19_time_series,. csv where confirm global is result and death is death global, just rename the two csv to respective .txt will do , with slight modification, the model can be run offline much more efficiently.

Mr Joshua, I am a Math's teacher and I want modeling the cases of Covid in my city Santa Marta, Colombia with his model. Can you tell me how i can do? thanks

Hi Joshua, please check on error for India, Indonesia, Japan, Singapore and Taiwan on examples.

Hi Joshu ,i have this error
fitVirusCV19v3("Italy",'prn','on')
Error using webread (line 122)
Could not establish a secure connection to "data.humdata.org". The reason is "error:14090086:SSL
routines:ssl3_get_server_certificate:certificate verify failed". Check your certificate file
(C:\Program Files\MATLAB\R2017b\sys\certificates\ca\rootcerts.pem) for expired, missing or invalid
certificates.

Joshua McGee

Will do.

St Se

Really thankful for this. Another examples update would be great.

Hi Joshu, could you please update data on examples tab regulary.

Ben Robes

Great, thanks! I've been admiring your model for quite some time!

Joshua McGee

This must be due to a typo in the HDX database. I'll look into this Ben.

Ben Robes

I believe your US cases per day for 4/21/20 is incorrect. From Wikipedia, the number is 25,304 (https://en.wikipedia.org/wiki/2020_coronavirus_pandemic_in_the_United_States). Your graph has the number close to 40,000.

Joshua McGee

The article is being updated.

Adit Negi

Can anyone please tell me how the getData funtion works? where is it pulling the data from? as I changed the data in the local csv files to implement it for indian states yet the results were the same i.e the results were from the old country data itself.

zhi li

Hello Mr.Joshua, this code is so great! But when I run it, I make frequent errors on the table index. I searched the help center and it said, "Starting with R2019b, you can specify table variable names that are not valid MATLAB® identifiers. Such variable names can include spaces, non-ASCII characters And there can be any character as the leading character. ", So only R2019b and above can be compatible with this program. Could you fix this compatibility?

Thanks so much Mr.Joshua. You've successfully modified Milan's SIR model. However, I can't see explicitly the equations of S, I and R. I searched in fitVirusCV19 & analyseCV19 but didn't find it. Can you let me know in which section you put these equations.

David Franco

Perfect Joshua! Thank you!

Joshua McGee

Hi David,

To fix this change ‘keeplimits’ to ‘keepticks’. This should do the trick.

Josh

Joshua McGee

Please keep in mind that until more data is available, the models will likely fail in providing accurate longterm predictions!

Jesse Fagan

I am getting similar error message (per Riel below) for Peru now:

fitVirusCV19v3("Peru")
Fail to obtain parameters.
ini: beta = 0.769831 gamma = 0.513221 N = 27777.4 I0 = 6.9504
calc: beta = 0.0292718 gamma = -0.169065 N = 602872 I0 = 112.561

Jesse Fagan

Joshua, that is great. Thank you so much. Here in Peru we are under a very strict quarantine since March 15th (not scheduled to end until Apr 26). I really appreciate your quick response.

Joshua McGee

Hi, the parameters beta and gamma are estimated off of the actual number of confirmed cases correct. Please remember that the variable R0 changes with respect to time. It will change every day based on the number of cases confirmed. While this number 2.4-2.5 is probably correct on average, it does not take into consideration the different control strategies that different countries are conducting. The model assumes equal contact between populations.
Josh

Jesse Fagan

Joshua, one more thing (and I may be getting way out of my league here!). Per the output for Peru, gamma = .7986 and beta = 1.286, R0 = 1.6103. Does the model estimate beta and gamma to fit the SIR model with the actual confirmed cases (number of infected)? It seems I am seeing a gamma rate much lower here in Peru (average .117). Also, for the coronavirus, in general, they are reporting a much higher average R0 value (2.0-2.5). Lastly, what assumptions does the model make about population mixing? Thanks again.

David Franco

Jesse Fagan

Many thanks, Joshua!

Joshua McGee

Hi Jesse,
Please see the example titled "R0s". This can be done by storing the output in an array from the command window. For example once in the same directory as the code type output = fitVirusCV19v3("US",'plt','off') and then output.R0 to see the R0 value and output.beta to see the beta value.
Josh

Jesse Fagan

Sorry, if this seems like s simple question. I noticed in the code that there are "optional outputs" for various parameters, but I can't seem to be able to get these numbers. What would the command look like for this output? I am interested in knowing res.RO, res.beta, res.gamma, and res.N. Thanks.

Thanks for the answer. Does the model take into account the total population for each country? If yes, where?

Joshua McGee

No, this is a fitted parameter. The max is a default setting included and can be changed. This is fitted as a part of the SIR model and therefore you don't need to touch this parameter. See attached: https://www.researchgate.net/publication/339311383_Estimation_of_the_final_size_of_the_coronavirus_epidemic_by_the_SIR_model

Why is Nmax = 2e6? Do I need to set it to the total population for each country?

Hi Joshua, please pay attention to the date report of the data. Example you mean, 1 day difference.

Joshua McGee

Hi Antonio,
I just investigated your claim. Here is what I found:
Predicted Total Cases:
Indonesia: (my code): 4387, (milan's code): 4446
Italy: (my code): 155889, (milan's code): 168747
US: (my code): 584464, (milan's code): 539470
Spain: (my code): 172017, (milan's code): 172105

Based on this I do not think that the difference that you are seeing is significant. This is definitely a result of differences how the data is being processed on both ends. It is expected that the values will be different and I do not believe that there is a huge significant difference. Thanks for looking into this.
Josh

Australia, Indonesia, Italy and many more. Try to check randomly.

Joshua McGee

Hi Antonn, which country are you looking at? Inherently there will be differences because the data sources are not the same.

Hi Joshua, based on same data as of April 5, there remains a significant difference with the results of Mr. Batista. Can you check it again?

Joshua McGee

Hi Gergely, I need more specifics than that. What is the error you are receiving?

Hello there. I'm trying to run your code, however, I'm getting "Error with country" - for all countries...

Thx Joshua, really really appreciate it...

Joshua McGee

Anton, I noticed there was a small mistake in the latest version of my data processing algorithm. Please try the most recent version. You should find the approximations to be more similar now.

Joshua McGee

I'll look into this

Joshua McGee

Hi Anton, the amount of data fed into the models may differ. Each code may select differing amounts of slow growth data and can be changed by simple modification of the code.

If this uses the same algorithm as Mr. Batista, why can the results be different (04/04 data)?

Aldrin Loly

Hi Mr. Joshua, I tried to run the code but it seems on getting an error of the code [t,Ca] = ode45(@(t,y) odeFun(t,y,b), tspan, I0) in line 338. Any idea why?

It doesn't work for Matlab 2010b.

yan binghai

Fantastic code!

Hi, Mr. Joshua! I succes run the it, but failed ( error ) to get some plot such as United Arab Emirates, Guam, New Caledonia.
On the otherhands the results for Turkey and Russia is good but the date axis for Turkey in case/day graph is unusual, the fit result of Russia show the large number of C_end.
Would you mind to explain these factors?
I am begging you.
Best regards!

vitto66

I am sorry, I just saw a similar question put on 30 march without reply

Roy CHEN

Hi, Mr. Joshua! It's very impresive work you done here!
I am from Hongkong( born in Taiwan). If I may ask, could you make update for my country Hongkong and Taiwan, please?
Your work will help us.
And please give the update for the exmaples (maybe) every week? ( If you can't )
Thank you very much for your kindness!

Thank you very much, Mr. Joshua McGee for updating my country Japan and the others!
Sorry to asking you a basic question. How can I install this model source?P.S. I used to use Python for my analytic.
I try to use it to see the other countries such as Russia, Vietnam, Turkey.
If you can tell me or update it for me, it would be great.
Thank you very much again! GBU

dong gretch

please give link of latest version

Joshua McGee

Hi Riel,
Thanks for reporting this bug. It should be fixed now. Please download the latest version.
Josh

Riel Aljama

Good morning Joshua,
I was able to run the your code now.
but I am having problem with some country, particularly Philippines

fitVirusCV19v3('Philippines','prn','on')
Fail to obtain parameters.
ini: beta = 0.682172 gamma = 0.454782 N = 8016.47 I0 = 0.00515606
calc: beta = 0.0543995 gamma = -0.219728 N = 25201.8 I0 = 0.000941286
Unrecognized function or variable 'plt'.

Error in fitVirusCV19v3>plotData (line 781)
if plt

Error in fitVirusCV19v3 (line 294)
plotData(C,date0,country)

Joshua McGee

Done!

Hi Mr Joshua McGee!
Thank you very much for all of these.
Could you do me a favor to add example case in Japan, my country?
Also if it is not trouble you, could you add Singapore and Australia too ?
Thank you very much. GBU

Joshua McGee

Hi Clarke! Thanks for your kind words. I'll work on incorporating per capita testing into my next iteration of the model.
Josh

Clarke

Hi Joshua, thank you very much for all this work. I'm wondering if you think increasing per capita COVID-19 testing rates here in the US (compared to other nations) will affect your total infected forecast?

Joshua McGee

Hi Marco, this code requires 2019b or greater.

Marco

Hi Joshua, this seems great. I am trying to run it but I get an error ... It seems like detectImportOptions is missing, so I am a bit stuck.
I am running on R2016a so it may be a bit dated.
Thanks for any tip
Marco

>> fitVirusCV19v3("Italy")
'lambertw' requires Symbolic Math Toolbox.

Error in fitVirusCV19v3>calcEndPoint (line 628)
lambertw(-beta*(1 - c0)*exp(-beta/gamma)/gamma);

Error in fitVirusCV19v3>calcClim (line 601)
res = calcEndPoint(beta,gamma,I0/N)*N;

Error in fitVirusCV19v3 (line 306)
Clim = calcClim(b);

Unfortunately, I don't have access to Symbolic Math Toolbox, and it's not listed as required to run the code. Any thoughts or suggestions?

Joshua McGee

added Indonesia to examples

Joshua McGee

Will do Anton.

Hi Joshua,
Would you mind making on examples tab for Indonesia? Thx.

Error using fitVirusCV19v3 (line 152)
Table variable names must be character vectors.
I'm using 2018 version of MATLAB..

Joshua McGee

-updated, attempted again to fix compatibility issues with older MATLAB versions

Joshua McGee

Will do. I’ll update later today.

St Se

Could you update the 'examples' graphs every couple of days maybe? Several countries approach the tipping point soon. Thanks

Joshua McGee

I’ll look into this.

on R2019b there is the following error:

at = fitVirusCV19v3("Austria");
Error using odearguments (line 21)
When the first argument to ode45 is a function handle, the tspan argument must have at least two elements.

Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in fitVirusCV19v3 (line 327)
[t,Ca] = ode45(@(t,y) odeFun(t,y,b), tspan, I0);

Thanks
Giovanni B.

Joshua McGee

It is fully functional for 2019b and up.

Joshua McGee

Yes this is a compatibility issue. I’m working on a workaround.

Hi! Cool model.
I keep geting:

Error using fitVirusCV19v3 (line 149)
'01_22_20' is not a valid variable name.

Error in run (line 1)
it = fitVirusCV19v3("Italy");

Any suggestions?

Joshua McGee

Hi Tosaporn, I’m working on a fix that will accommodate most MATLAB users. This is compatible with 2019b and up.

Hi Joshua,
I tried the code, it showed Error , Can you please explain?
at = fitVirusCV19v3("Austria");
Error using fitVirusCV19v3 (line 152)
'01_22_20' is not a valid table variable name.
See the documentation for isvarname or
matlab.lang.makeValidName for more information.

Joshua McGee

Hi Giuseppe,
The database was probably experiencing too much traffic at the time. Wait a little while and try again.
Josh

Leon Burger

Hi, these are my results on R2019a; any work around. Want to run scenario for South Africa, but using Italty as a test case.

Error using readContentFromWebService (line 46)
The server returned the status 403 with message "FORBIDDEN" in response to the request to URL
https://proxy.hxlstandard.org/api/data-preview.csv?url=https%3A%2F%2Fraw.githubusercontent.com%2FCSSEGISandData%2FCOVID-19%2Fmaster%2Fcsse_covid_19_data%2Fcsse_covid_19_time_series%2Ftime_series_19-covid-Confirmed.csv&filename=time_series_2019-ncov-Confirmed.csv&options=table.

Error in webread (line 125)
[varargout{1:nargout}] = readContentFromWebService(connection, options);

Error in fitVirusCV19v3 (line 124)
result=webread('https://proxy.hxlstandard.org/api/data-preview.csv?url=https%3A%2F%2Fraw.githubusercontent.com%2FCSSEGISandData%2FCOVID-19%2Fmaster%2Fcsse_covid_19_data%2Fcsse_covid_19_time_series%2Ftime_series_19-covid-Confirmed.csv&filename=time_series_2019-ncov-Confirmed.csv','options','table');

Joshua McGee

Hi Riel,
Thanks for your patience. Please try again now. I converted the variable name dates to datenum format.
Josh

Riel Aljama

Hi, these are the result

Error using fitVirusCV19v3 (line 148)
'22-Jan-2020' is not a valid variable name.

Joshua McGee

Hi please try again now.

Riel Aljama

Good Day!
I keep on getting this error :( (using 2018a)

Error using fitVirusCV19v3 (line 148)
'Province/State' is not a valid variable name.

Joshua McGee

All compatibility issues should now be fixed

Joshua McGee

This is a compatibility issue. What version are you running?

any ideas on why I would get this error:

Unrecognized property 'PreserveVariableNames' for class
'matlab.io.text.DelimitedTextImportOptions'.

Error in fitVirusCV19v3 (line 128)
opts.PreserveVariableNames = true;

St Se

Very nice to see the graphs being updated.

St Se

Suksan,
I had the same problem. Matlab 2019.b with mex installed.

excellent work, but when I run the code it gives me this error:
Unrecognized property 'PreserveVariableNames' for class
'matlab.io.text.DelimitedTextImportOptions'.
Error in fitVirusCV19v3 (line 132)
opts.PreserveVariableNames = true;

buldus2

excellent work, but when I run the code it gives me this error: Error using fitVirusCV19v3 (line 170)
Duplicate table variable name: 'NaN'.

any suggestion?

Updates

3.4.8

updated 5/2

3.4.7

updated 5/1/20

3.4.6

updated 4/30

3.4.5

updated 4/22

3.4.4

updated 4/22

3.4.3

updated info

3.4.2

updated 4/18

3.4.1

updated 4/17

3.4

-updated notes and data for 4/15

3.3.9

updated with link to in-depth description of model

3.3.8

updated 4/15

3.3.7

updated 4/13

3.3.6

-updated

3.3.5

-slight bug fixes (error with plotting data that failed to fit the model)
-adjusted accuracy by changing data selection algorithm

3.3.4

updated 4/12, introduced the ability to plot changes in R0 over time with 'r0analysis' command

3.3.3

updated 4/11

3.3.2

-updated 4/10

3.3.1

updated 4/9

3.3

-updated 4/8

3.2.9

-updated (4/6)

3.2.8

-updated

3.2.7

-small fix

3.2.6

-image

3.2.5

-minor updates and corrections

3.2.4

-improved initial parameter approximation (courtesy of Milan Batista)

3.2.3

-updated examples

3.2.2

-updated with data from 4/2

3.2.1

-image

3.2

-updated 4/1

3.1.9

-bug fixes

3.1.8

-updated 3/31 8:40 PM EST

3.1.7

-updated aesthetics of examples

3.1.6

-updated examples

3.1.5

-updated plots for today (3/30)

3.1.4

-bug fix

3.1.3

-updated examples

3.1.2

-image

3.1.1

-compatibility notes

3.1

-updated examples

3.0

-updated plots and fixed small bug

2.9

-added death rate to example

2.8

-updated examples

2.7

-added ability to disable plot and created a live script to tabulate the R0 values of each country

2.6

-added Indonesia

2.5

-bug fixes and compatibility

2.3

-compatibility

2.2

-attempted to fix additional variable name Compatability issues

2.1

-compatibility

2.0.9

-fixed all compatibility issues

2.0.8

-bug fixes

2.0.7

-image

2.0.6

-description

2.0.5

-updated date ticks

2.0.4

-new data source (updated database)

2.0.3

-description

2.0.2

-updated compatibility

2.0.1

-updated description

2.0

-applied updates courtesy of Milan Batista (RMSE correction, automatic parameter selection)

1.0.4

-name change

1.0.3

-added tags

1.0.2

-fixed live script

1.0.1

-added live script to generate plots

MATLAB Release Compatibility
Created with R2019b
Compatible with R2019b to any release
Platform Compatibility
Windows macOS Linux