alsnhll / seir_covid19 Goto Github PK
View Code? Open in Web Editor NEWSEIR model for COVID-19 infection, including different clinical trajectories of infection
SEIR model for COVID-19 infection, including different clinical trajectories of infection
Hi @alsnhll ,
I am investigating the accuracy of the results following an intervention as implemented in the app.
My test is the following:
I want to run an intervention with zero effect. The implementation is accurate if the projections following the intervention are identical with the results in the base scenario (baseline). We can do this by setting the s0, s1, s2, s3
parameters equal to zero.
Find below a reproducible example mostly using the code in runIntervention.R
Define the inputs and call the functions:
> remove(list = ls())
>
> library(deSolve)
> library(plotly)
> library(dplyr)
> library(reshape)
> library(htmltools)
> source("code/functions.R")
>
> # Set parameters
> # (in Shiny app these are input by sliders)
>
> IncubPeriod= 5 #Duration of incubation period
> DurMildInf= 6 #Duration of mild infections
> FracSevere= 15 #% of symptomatic infections that are severe
> FracCritical= 5 #% of symptomatic infections that are critical
> ProbDeath= 40 #Death rate for critical infections
> DurHosp= 5 #Duration of severe infection/hospitalization
> TimeICUDeath= 8 #Duration critical infection/ICU stay
> AllowPresym="No"
> AllowAsym="No"
> FracAsym=25 #Fraction of all infections that are asymptomatic
> PresymPeriod=2 #Length of infectious phase of incubation period
> DurAsym=6 #Duration of asympatomatic infection
> b1= 0.5 #Transmission rate (mild infections)
> b2= 0.1 #Transmission rate (severe infections)
> b3= 0.1 #Transmission rate (critical infections)
> be = 0.5 #Transmission rate (pre-symptomatic)
> b0 = 0.5 #Transmission rate (asymptomatic infections)
> N= 1000 #Total population size
> Tmax= 300 #Initial # infected
> InitInf= 1 #Maximum time
> yscale="linear"
> AllowSeason="No"
> seas.amp=0.0 #relative amplitude of seasonal fluctuations, in [0,1]
> seas.phase=0 #phase of seasonal fluctuations, measuered in days relative to time zero when peak will occur (0=peak occurs at time zero, 30 = peak occurs one month after time zero). Can be negative
> #Note that the beta values chosen above will be those for time zero, wherever that occurs in the seasonality of things.
>
> #Intervention parameters
> VarShowInt="E" #This is the variable to be plotted. Options "E0", "E1","I0", "I1", "I2", "I3", "R", "D", "Suceptible (S)"="S", "Exposed (E)"="E", "Mild Infections (I1)"="I1", "Severe Infections (I2)"="I2", "Critical Infections (I3)"="I3", "Recovered (R)"="R", "Dead (D)"="D", "All infected (E + all I)"="Inf","All symptomatic (I1+I2+I3)"="Cases","All hospitalized (I2+I3)"="Hosp"),
> Tint = 30 #Intervention start time (days)
> Tend = 300 #Intervention end time (days)
>
> # NOTE ALL THE REDUCTION PARAMETERS ARE EQUAL TO ZERO
> # THIS MEANS AN INTERVENTION WITH ZERO EFFECT
> #
> s0 = 0 #Reduction in transmission from asymptomatic/presymptomatic infections
> s1 = 0 #Reduction in transmission from mild infections
> s2 = 0 #Reduction in transmission from severe infections
> s3 = 0 #Reduction in transmission rate from critical infections
> RoundOne = FALSE #Round values to nearest integar post-intervention?"
>
> #Put these into an input structure
> input=list("IncubPeriod"=IncubPeriod,"DurMildInf"=DurMildInf,"FracSevere"=FracSevere,"FracCritical"=FracCritical,"ProbDeath"=ProbDeath,"DurHosp"=DurHosp,"TimeICUDeath"=TimeICUDeath,"FracAsym"=FracAsym, "PresymPeriod"=PresymPeriod, "DurAsym"=DurAsym, "be"=be,"b0"=b0,"b1"=b1,"b2"=b2,"b3"=b3,"seas.amp"=seas.amp, "seas.phase"=seas.phase,"N"=N,"Tmax"=Tmax,"InitInf"=InitInf,"yscale"=yscale,"AllowPresym"=AllowPresym,"AllowAsym"=AllowAsym,"AllowSeason"=AllowSeason,"VarShowInt"=VarShowInt,"Tint"=Tint,"Tend"=Tend,"s0"=s0,"s1"=s1,"s2"=s2,"s3"=s3,"RoundOne"=RoundOne)
The test here:
# Run simulations
sim = SimSEIR(input)
simInt = SimSEIRintB(input)
Check the dimensions of the output. We expect the same output, but it does not seem to be the case. But we already know this from #11. This should be an easy fix.
> dim(sim$out.df)
[1] 301 10
> dim(simInt$out.df)
[1] 302 10
Now, even if we remove the doubled row we still see differences:
> testthat::expect_identical(sim$out.df, simInt$out.df[-31, ])
Error: sim$out.df not identical to simInt$out.df[-31, ].
Attributes: < Component “row.names”: Mean relative difference: 0.006024096 >
Component “S”: Mean relative difference: 7.208045e-07
Component “E0”: Mean relative difference: 1.464121e-06
Component “E1”: Mean relative difference: 1.464117e-06
Component “I1”: Mean relative difference: 1.475496e-06
Component “I2”: Mean relative difference: 1.526325e-06
Component “I3”: Mean relative difference: 1.385202e-06
Component “R”: Mean relative difference: 1.056044e-07
Component “D”: Mean relative difference: 1.179242e-07
And this looks like this in the output data frame:
> (sim$out.df - simInt$out.df[-32, ])[29:34, ]
time S E0 E1 I0 I1 I2 I3 R D
29 0 0.000000e+00 0.000000e+00 0.000000e+00 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
30 0 0.000000e+00 0.000000e+00 0.000000e+00 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
31 0 0.000000e+00 0.000000e+00 0.000000e+00 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
32 0 5.455557e-05 -2.489039e-05 -4.973590e-12 0 -1.591729e-05 -1.388994e-06 -2.123629e-07 -1.209542e-05 -5.111330e-08
33 0 1.249899e-04 -4.452048e-05 -8.905145e-12 0 -3.520897e-05 -3.747634e-06 -7.147822e-07 -4.055353e-05 -2.445333e-07
34 0 1.721000e-04 -7.518219e-05 -1.503644e-11 0 -4.866380e-05 -4.605636e-06 -7.638421e-07 -4.268553e-05 -1.989793e-07
Let me know if I am missing something. Thanks a lot!
Marius
Great job, Alison! Thanks a lot for your work.
I am exploring the possibility of modifying the app to allow for multiple interventions of different intensity (like what we see in real world) and also to export the result in a Excel file.
Checking the output of the SimSEIRintB()
I can see that the day of the intervention is displayed twice (see below). Same values and I guess plotted like this.
Replicate the results with the code in runIntervention.R
> Tint = 30 # Intervention start time (days)
> simInt = SimSEIRintB(input)
>
> outInt.df = simInt$out.df
>
> outInt.df[30:33, ]
time S E0 E1 I0 I1 I2 I3 R D
30 29 942.5220 22.98966 4.597931e-06 0 15.53003 1.560992 0.3010213 16.99366 0.1026588
31 30 934.5564 26.05645 5.211289e-06 0 17.66609 1.779912 0.3441635 19.47826 0.1187653
32 30 934.5564 26.05645 5.211289e-06 0 17.66609 1.779912 0.3441635 19.47826 0.1187653
33 31 928.2269 27.08032 5.416064e-06 0 19.84364 2.025373 0.3931783 22.29337 0.1371733
Marius
Would it be possible to save or input all parameters in a JSON file? It's difficult to enter all values all the time over an over again. Thanks for this great tool!
With the current version of googlesheets4 (1.0.0), the sheets_deauth
function doesn't exist. The app will run if you define
sheets_deauth<-googlesheets4::gs4_deauth
We would like to use the John Hopkins data to estimate early stage exponential rate r. Do we have criteria for which time subset we would like to use for that? We could use something like the first ten days after fifty cases are reported, or something similarly arbitrary?
Shouldn't R0 be independent of N ?
If I have ten times the population, I should not have ten times the R0, correct?
Or, ... what am I missing?
This comment is based on the code in the Python notebook.
EDIT:
thanks.
Dear @alsnhll,
Thanks you for shared your work.
I am trying to understand why when we increase Reduction of transmission percentage, the impact of the intervention disappear and the two curves seems to be same with just a lag in time?
This occurs when the Intervention curve exceeds the Tend
.
With default setting I just reduce Tend
to 100 instead 300. And tried to increase Reduction in transmission from mild infections
to 70%.
Which line of code I have to understand?
Thanks
Thanks for making this great app available! The intervention and capacity tabs currently crash when adjusting the maximum time below the default 300 without also adjusting the intervention end time accordingly.
Capping the intervention end time at the maximum time of the simulation solves this. Or just replacing
Tend=input$Tend
with
Tend=pmin(input$Tend,input$Tmax)
in the SimSEIRintB
function.
If the Initial # infected:
slider is set tot 1000, then the Capacity tab produces this result.
Reading from 'COVID19 Parameters'
Range "'Hospital capacity'"
New names:
* `` -> ...6
Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(V)` instead of `V` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.
Warning in max(abs(diff(times))) :
no non-missing arguments to max; returning -Inf
Warning: Error in checkInput: `hmax' must be a non-negative value
104: stop
103: checkInput
102: lsoda
101: ode
100: GetSpread_SEIR [/media/truecrypt1/Projects/SEIR_COVID19/COVID19seir/server.R#49]
99: renderPlotly [/media/truecrypt1/Projects/SEIR_COVID19/COVID19seir/server.R#493]
98: func
95: func
82: origRenderFunc
81: output$plotCap
1: runApp
When I try to increase the N from 1000 to 100,000; The plot changes dramatically, which it should not. Did I miss anything ? Thanks
Hey,
So while running the capacity file, there is a error I get, as follows:
Run simulations
sim=SimSEIR(input)
Error in if (input$AllowSeason == "Yes") { : argument is of length zero
Can someone help with the resolution of the same
Thank you.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.