Intervention procedure issue

[mpascariu]

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

[mpascariu]

Maybe the errors are negligible. I was afraid of a snowball effect.