alsnhll / seir_covid19

SEIR model for COVID-19 infection, including different clinical trajectories of infection

SEIR_COVID19

The code for creating the R Shiny application https://alhill.shinyapps.io/COVID19seir/ is in the directory COVID19seir. This code is an SEIR model for COVID-19 infection, including different clinical trajectories of infection, interventions to reduce transmission, and comparisons to healthcare capacity.

The code that produces the interface and functionality of the Shiny App is in files

• server.R
• ui.R

Files used in the explanatory sections of the app are

• SEIR.Rmd
• About.Rmd
• www/Parameters.nb.html
• www/Tutorial.html

All the functions that actually run the model and process the parameters are in the code/functions.R file

If you want to run the code to produce the same outputs as Shiny but without dealing with the app structure, you can use the R scripts

• runSpread.R
• runInterventions.R
• runCapacity.R

When trying out new model structures or plots, it is much easier to work with scripts instead of directly with the app. Once troubleshooting is done, it can be integrated with the app

There are also some new new scripts to simulate the worsening of outcomes if healthcare capacity is overwhelmed

• runOverflow.R
• code/functionsOverflow.R

There is also a Python notebook simulating the same model (.ipynb)

You can use this page to report issues you find with the app

Developers: No pull requests will be accepted from this repository, which is only for sharing the stable code used in the app. If you are contributing to this project please use the developmental repository https://github.com/alsnhll/SEIR_COVID19_Dev

Caution: This model is still a work in progress and will be updated continuously.

seir_covid19's Issues

Save/Input all parameters at once

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!

Error: `hmax' must be a non-negative value

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

ODE chart correlated to N

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

Estimation of r

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?

Intervention day saved twice

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

Intervention procedure issue

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

New version of googlesheets4

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

why is R, the rate of spread, calculated as a multiplier of N, the total population size?

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:

• changed R to R0 in the description and text, realizing I was initially confusing one for the other.
• deleted initial context, unnecessary

thanks.

Intervention plots return errors at maximum times <100 days

Intervention when Reduction in transmission from mild infections is bigger that 30%

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

EpiModel issue

Problem with maximum time and intervention end time mismatch

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.

Problem with Capacity File

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.