Nothing
Stochastic SEIR Model
In DSAIDE: Dynamical Systems Approach to Infectious Disease Epidemiology (Ecology/Evolution)
#load various variable definitions that are the same for each app
source('startup_script.R')
currentrmdfile = knitr::current_input()
appsettings = get_settings(currentrmdfile,appdocdir,packagename)
Overview {#shinytab1}
This app lets you explore a stochastic SEIR model. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.
For this app, it is assumed that you've worked through all the ones in the Basics and Reproductive Number sections, as well as the Stochastic SIRS one.
Learning Objectives
- Obtain further familiarity with stochastic compartmental models
- Learn about steady states in stochastic models
- Explore extinctions
The Model {#shinytab2}
Model Overview
This model tracks susceptibles, exposed/pre-symptomatic, infected/symptomatic and recovered hosts.
The following compartments are included:
- S - uninfected and susceptible individuals.
- E - individuals who are infected and do not yet show symptoms. Those individuals can be infectious.
- I - individuals who are infected and show symptoms. Those individuals are likely infectious, but the model allows to adjust this, including no infectiousness.
- R - recovered/removed individuals. Those individuals have recovered and are immune. They can lose their immunity and return to the susceptible class.
The included processes/mechanisms are the following:
- Susceptible individuals (S) can become infected by pre-symptomatic/exposed (E) or symptomatic (I) hosts at rates b~E~ and b~I~.
- All infected hosts first enter the exposed/presymptomatic stage and leave the stage at rate g~E~, the inverse of which is the average time spent in this stage.
- All infected hosts go from the presymptomatic to the symptomatic stage and leave the stage at rate g~I~.
- Recovered individuals are initially immune to reinfection. They can lose their immunity at rate w and return to the susceptible compartment.
- New susceptible hosts enter the system (are born) at fixed rate n. Natural death occurs at rate m.
Model Implementation
The flow diagram for the model implemented in this app is:
knitr::include_graphics(here::here('inst/media',appsettings$modelfigname))
The deterministic/ODE equations for this model are:
$$\dot S = n - S (b_E E + b_I I) - mS + wR$$
$$\dot E = S (b_E E + b_I I) - g_E E - m E$$
$$\dot I = g_E E - g_I I - m I $$
$$\dot R = g_I I - mR - wR$$
However, for this app we do not implement a deterministic/ODE model. Instead, we implement its stochastic equivalent.
We can specify the model by writing down every possible transition/event/reaction that can occur and their propensities (the propensity multiplied with the
time step gives the probability that a given event/transition occurs). For our model these are the following:
Event type | Transitions | Propensity |
---------- | ----------- | ---------- |
Infection | S => S-1, E => E+1 | S(b~E~E+b~I~I) |
Progression to Symptoms | E => E-1, I => I+1 | g~E~E |
Recovery | I => I-1, R => R+1 | g~I~I |
Waning of Immunity | R => R-1, S => S+1 | wR |
Births | S => S+1 | n |
Death of susceptible | S => S-1 | mS |
Death of exposed | E => E-1 | mE |
Death of symptomatic | I => I-1 | mI |
Death of recovered | R => R-1 | mR |
What to do {#shinytab3}
The tasks below are described in a way that assumes everything is in units of MONTHS (rate parameters, therefore, have units of inverse months). If any quantity is not given in those units, you need to convert it first (e.g. if it says a year, you need to convert it to 12 months).
Some of the simulations might take a few seconds to run. Be patient.
#this is the running counter for the records which starts at 1
rc=1
#empty object, will hold all outcomes
alloutcomes = NULL
#########################
# Task 1
#########################
tid = 1
tasktext = "Set the following model parameters to start 1000 susceptible and 1 initially symptomatic host, none in the other categories. Assume that only symptomatic individuals transmit, at rate 0.005. Assume that the duration of the presymptomatic period is one week (1/4 of a month) long, and the duration of the symptomatic period is 2 weeks (half a month) long. Assume immunity does not wane and that there are no births and deaths. Set the number of simulations to 10. Simulation start time 0, duration 1 year, time step does not matter. Set random seed to 123. With parameters set to correspond to the scenario just described, run the simulation. You should find 5 simulations that produced an outbreak and the average final value for the recovered is around 448."
nrec = 1 # number of items to record
out_records = c("Average number of susceptible at end of simulation")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#########################
# Task 2
#########################
tid = tid + 1
tasktext = "Now set waning immunity rate to 0.2, simulation time 5 years. Keep everything as before. Run the simulation, you should find that the system settles to some kind of stochastic version of a steady state, namely S-E-I-R values that fluctuate around some level. Note that there is still a mix of simulations that did and didn't produce an outbreak. In fact, that number of simulations that produce an outbreak dropped from 5 to 4. Just by chance the inclusion of waning immunity meant that processes in the model got executed differently at the start of the simulation and a simulation that took off earlier now didn't. Let's - just for fun - look into that more closely and see if we can identify the simulation. To that end, set the number of simulations to run to 1, start with random seed 123, no waning immunity. Run the simulation. You'll see that this produces no outbreak. We want to find the one that produces an outbreak in the absence of waning immunity and does not with waning immunity turned on. So let's move to the next one, run the simulation for random seed 124 (remember from the previous app that if you run 10 simulations, internally the seed is increased each time). Again no outbreak, so let's do 125. You'll see an outbreak. Now turn on waning immunity see if you still get an outbreak (you will). Keep trying the random seeds (up to 132) until you find the seed for which you get an outbreak without waning immunity but none with waning immunity. Note that this has nothing directly to do with the waning immunity process, it just happens to mix up the random numbers in the simulation that now we don't get an outbreak."
nrec = 1 # number of items to record
out_records = c("Random number seed for which you get an outbreak in the absence of waning immunity and no outbreak in the presence of waning immunity")
out_types = rep("Integer",nrec)
out_notes = rep("Report the integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#########################
# Task 3
#########################
tid = tid + 1
tasktext = "Let's go back to the settings at the beginning of the previous task, namely w=0.2, 10 simulations, random seed 123. You can see that the reported final value for the susceptibles does not agree with what you see for those simulations which have a steady-state like scenario with non-zero infected. That's because the average includes also those scenarios where no outbreak happened. Essentially, the average is computed here by doing 6/10\\*1000 (simulations without an outbreak) + 4/10\\*Sfinal_outbreak (simulations with outbreak), where Sfinal_outbreak is the average value of susceptibles at the end of the simulation for those scenarios where you get an outbreak. Approximately read off a value for Sfinal_outbreak and convince yourself that adding the non-outbreak and outbreak scenarios together gives you the reported mean."
nrec = 1 # number of items to record
out_records = c("Average number (across all simulations) of susceptible at end of simulation, as reported by DSAIDE")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#########################
# Task 4
#########################
tid = tid + 1
tasktext = "Sometimes you do want to know the overall average, including both settings where an outbreak happened and where it didn't. However, often combining those two very different scenarios can be confusing. So let's change our setup to make sure we get outbreaks every time. You learned previously that if there are more initial infected, chances of an outbreak are larger. Let's do that. Keep all settings as before, but change initial number of symptomatic infected to 2. You'll see that there are a few more outbreaks than with 1 infected, but still not all. Keep increasing the starting value for _I_ until all 10 simulations lead to outbreaks."
nrec = 1 # number of items to record
out_records = c("Minimum number infected such that all 10 simulations produce outbreaks")
out_types = rep("Integer",nrec)
out_notes = rep("Report the integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#########################
# Task 5
#########################
tid = tid + 1
tasktext = "Set initial number infected to 10, set gI=1, leave all other settings as in the previous task (bE=0, bI=0.005, gE=4, w=0.2). Run the simulation for 5 years. You should find an average of around 133 infected at the end. While there is no guarantee that the average of the stochastic model agrees with an equivalent deterministic model for the types of (nonlinear) models we are exploring, they are often close. We can check that by computing the steady states for the deterministic version of the model. If you want to practice doing steady states, you can go ahead and solve the SEIR ODE model for the steady states of the variables. See e.g. the _Patterns of ID_ app for an explanation on how to do it. If you don't want to practice, I'll tell you that the steady state value for S is S=(g~E~ * g~I~)/(b~E~ * g~I~ + b~I~ * g~E~). Use the current values to compute S and compare with results from the simulation."
nrec = 2 # number of items to record
out_records = c("Theoretical value for susceptible at steady state, based on equation",
"Average number (across all simulations) of susceptible at end of simulation")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#########################
# Task 6
#########################
tid = tid + 1
tasktext = "Let's repeat the previous simulation, but turn off waning immunity. Also set back to gI=2. Still 10 initial infected, random seed at 123, do 20 simulations (be patient), run for 5 years. Record the average number of susceptible left at the end of the simulation. Use that value to compute R~0~ using the final size equation. Then compare it with the R~0~ value computed based on the model (this is the same as for the SIR model, the E compartment does not make a difference here)."
nrec = 2 # number of items to record
out_records = c("Average number (across all 20 simulations) of susceptible at end of simulation",
"R0 based on model equation")
out_types = c("Rounded_Integer","Numeric")
out_notes = c("Report the rounded integer","Report to one decimal place")
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#########################
# Task 7
#########################
tid = tid + 1
tasktext = "Finally, let's explore extinctions. We'll start with 100 individuals in each compartment, b~E~=0 and b~I~=0.01, g~E~=4, g~I~=2, w=1. No births and deaths, start simulation at 0, run for 60 months, time step arbitrary. Run 10 simulations with seed of 123. You should see some initial ups and downs, and then for each simulation, you'll see fluctuations around some type of steady state, with an average number of infected at the end of 57. For those settings, in the absence of transmission from the _E_ class, the basic reproductive number (if everyone were susceptible) is R~0~ = S~0~b~I~/g~I~ = 0.01*400/2 = 2. Now, reduce transmission to bI=0.006. Compute R~0~ for this value. Run simulations for 60 months. You will find that for some of these simulations, the pathogen goes extinct (and susceptible go up to the maximum value of 400). Now double simulation time to 10 years. You will find that for most simulations, extinction occured."
nrec = 2 # number of items to record
out_records = c("Number of simulations that DO NOT lead to extinction after 10 years of simulation (with bI=0.006)",
"Reproductive number for bI=0.006")
out_types = c("Rounded_Integer","Numeric")
out_notes = c("Report the rounded integer","Report to one decimal place")
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#########################
# Task 8
#########################
tid = tid + 1
tasktext = "Keep everything as in the previous task, but change infected tranmission rate to b~I~=0.007. Compute R~0~ for this value. Run 10 simulations for 10 years. You will find that for most simulations, extinction did not occur."
nrec = 2 # number of items to record
out_records = c("Number of simulations that DO lead to extinction after 10 years of simulation (with bI=0.007)",
"Reproductive number for bI=0.007")
out_types = c("Rounded_Integer","Numeric")
out_notes = c("Report the rounded integer","Report to one decimal place")
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#########################
# Task 9
#########################
tid = tid + 1
tasktext = "Keep exploring. You can test how R~0~ and other model quantities affect extinction, and how the longer you run the simulations, the more extinctions you get. You can also turn on births and deaths and see how they affect the overall patterns you find. There is a lot more to explore. For some settings, i.e. if you want to run many replicates for a long time, things are starting to slow down. That's a disadvantage of stochastic models, they often take longer to run."
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
TaskText = rep(tasktext,nrec),
RecordID = paste0('T',tid,'R',(1:nrec)),
Record = out_records,
Type = out_types,
Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task
#save the fully filled task table to a tsv file
alloutcomes$QuizID = paste0(packagename,"_",appsettings$appid)
alloutcomes$AppTitle = appsettings$apptitle
alloutcomes$AppID = appsettings$appid
#remove a few variables from the data frame
savedoutcomes <- dplyr::select(alloutcomes,QuizID,AppID,AppTitle,TaskID,TaskText,RecordID,Record,Type,Note)
write.table(savedoutcomes, paste0(appsettings$appid,"_tasktable.tsv"), append = FALSE, sep = "\t", row.names = F, col.names = TRUE)
# Take all the text stored in the table and print the tasks and items to record
write_tasktext(alloutcomes)
Further Information {#shinytab4}
- This app (and all others) are structured such that the Shiny part (the graphical interface you see and the server-side function that goes with it) calls an underlying R script (or several) which runs the simulation for the model of interest and returns the results.
- For this app, the underlying function(s) running the simulation are called
r appsettings$simfunction. You can call them directly, without going through the shiny app. Use the help() command for more information on how to use the functions directly. If you go that route, you need to use the results returned from this function and produce useful output (such as a plot) yourself.
- You can also download all simulator functions and modify them for your own purposes. Of course to modify these functions, you'll need to do some coding.
- For examples on using the simulators directly and how to modify them, read the package vignette by typing
vignette('DSAIDE') into the R console.
- For some further reading concerning stochastic dynamics of ID, see e.g. [@lloyd-smith05; @black1966; @keeling1997].
References
Try the DSAIDE package in your browser
Any scripts or data that you put into this service are public.
DSAIDE documentation built on Aug. 24, 2023, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#load various variable definitions that are the same for each app source('startup_script.R') currentrmdfile = knitr::current_input() appsettings = get_settings(currentrmdfile,appdocdir,packagename)
Overview {#shinytab1}
This app lets you explore a stochastic SEIR model. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.
For this app, it is assumed that you've worked through all the ones in the Basics and Reproductive Number sections, as well as the Stochastic SIRS one.
Learning Objectives
- Obtain further familiarity with stochastic compartmental models
- Learn about steady states in stochastic models
- Explore extinctions
The Model {#shinytab2}
Model Overview
This model tracks susceptibles, exposed/pre-symptomatic, infected/symptomatic and recovered hosts. The following compartments are included:
- S - uninfected and susceptible individuals.
- E - individuals who are infected and do not yet show symptoms. Those individuals can be infectious.
- I - individuals who are infected and show symptoms. Those individuals are likely infectious, but the model allows to adjust this, including no infectiousness.
- R - recovered/removed individuals. Those individuals have recovered and are immune. They can lose their immunity and return to the susceptible class.
The included processes/mechanisms are the following:
- Susceptible individuals (S) can become infected by pre-symptomatic/exposed (E) or symptomatic (I) hosts at rates b~E~ and b~I~.
- All infected hosts first enter the exposed/presymptomatic stage and leave the stage at rate g~E~, the inverse of which is the average time spent in this stage.
- All infected hosts go from the presymptomatic to the symptomatic stage and leave the stage at rate g~I~.
- Recovered individuals are initially immune to reinfection. They can lose their immunity at rate w and return to the susceptible compartment.
- New susceptible hosts enter the system (are born) at fixed rate n. Natural death occurs at rate m.
Model Implementation
The flow diagram for the model implemented in this app is:
knitr::include_graphics(here::here('inst/media',appsettings$modelfigname))
The deterministic/ODE equations for this model are:
$$\dot S = n - S (b_E E + b_I I) - mS + wR$$ $$\dot E = S (b_E E + b_I I) - g_E E - m E$$ $$\dot I = g_E E - g_I I - m I $$ $$\dot R = g_I I - mR - wR$$
However, for this app we do not implement a deterministic/ODE model. Instead, we implement its stochastic equivalent. We can specify the model by writing down every possible transition/event/reaction that can occur and their propensities (the propensity multiplied with the time step gives the probability that a given event/transition occurs). For our model these are the following:
Event type | Transitions | Propensity | ---------- | ----------- | ---------- | Infection | S => S-1, E => E+1 | S(b~E~E+b~I~I) | Progression to Symptoms | E => E-1, I => I+1 | g~E~E | Recovery | I => I-1, R => R+1 | g~I~I | Waning of Immunity | R => R-1, S => S+1 | wR | Births | S => S+1 | n | Death of susceptible | S => S-1 | mS | Death of exposed | E => E-1 | mE | Death of symptomatic | I => I-1 | mI | Death of recovered | R => R-1 | mR |
What to do {#shinytab3}
The tasks below are described in a way that assumes everything is in units of MONTHS (rate parameters, therefore, have units of inverse months). If any quantity is not given in those units, you need to convert it first (e.g. if it says a year, you need to convert it to 12 months).
Some of the simulations might take a few seconds to run. Be patient.
#this is the running counter for the records which starts at 1 rc=1 #empty object, will hold all outcomes alloutcomes = NULL ######################### # Task 1 ######################### tid = 1 tasktext = "Set the following model parameters to start 1000 susceptible and 1 initially symptomatic host, none in the other categories. Assume that only symptomatic individuals transmit, at rate 0.005. Assume that the duration of the presymptomatic period is one week (1/4 of a month) long, and the duration of the symptomatic period is 2 weeks (half a month) long. Assume immunity does not wane and that there are no births and deaths. Set the number of simulations to 10. Simulation start time 0, duration 1 year, time step does not matter. Set random seed to 123. With parameters set to correspond to the scenario just described, run the simulation. You should find 5 simulations that produced an outbreak and the average final value for the recovered is around 448." nrec = 1 # number of items to record out_records = c("Average number of susceptible at end of simulation") out_types = rep("Rounded_Integer",nrec) out_notes = rep("Report the rounded integer",nrec) outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task ######################### # Task 2 ######################### tid = tid + 1 tasktext = "Now set waning immunity rate to 0.2, simulation time 5 years. Keep everything as before. Run the simulation, you should find that the system settles to some kind of stochastic version of a steady state, namely S-E-I-R values that fluctuate around some level. Note that there is still a mix of simulations that did and didn't produce an outbreak. In fact, that number of simulations that produce an outbreak dropped from 5 to 4. Just by chance the inclusion of waning immunity meant that processes in the model got executed differently at the start of the simulation and a simulation that took off earlier now didn't. Let's - just for fun - look into that more closely and see if we can identify the simulation. To that end, set the number of simulations to run to 1, start with random seed 123, no waning immunity. Run the simulation. You'll see that this produces no outbreak. We want to find the one that produces an outbreak in the absence of waning immunity and does not with waning immunity turned on. So let's move to the next one, run the simulation for random seed 124 (remember from the previous app that if you run 10 simulations, internally the seed is increased each time). Again no outbreak, so let's do 125. You'll see an outbreak. Now turn on waning immunity see if you still get an outbreak (you will). Keep trying the random seeds (up to 132) until you find the seed for which you get an outbreak without waning immunity but none with waning immunity. Note that this has nothing directly to do with the waning immunity process, it just happens to mix up the random numbers in the simulation that now we don't get an outbreak." nrec = 1 # number of items to record out_records = c("Random number seed for which you get an outbreak in the absence of waning immunity and no outbreak in the presence of waning immunity") out_types = rep("Integer",nrec) out_notes = rep("Report the integer",nrec) outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task ######################### # Task 3 ######################### tid = tid + 1 tasktext = "Let's go back to the settings at the beginning of the previous task, namely w=0.2, 10 simulations, random seed 123. You can see that the reported final value for the susceptibles does not agree with what you see for those simulations which have a steady-state like scenario with non-zero infected. That's because the average includes also those scenarios where no outbreak happened. Essentially, the average is computed here by doing 6/10\\*1000 (simulations without an outbreak) + 4/10\\*Sfinal_outbreak (simulations with outbreak), where Sfinal_outbreak is the average value of susceptibles at the end of the simulation for those scenarios where you get an outbreak. Approximately read off a value for Sfinal_outbreak and convince yourself that adding the non-outbreak and outbreak scenarios together gives you the reported mean." nrec = 1 # number of items to record out_records = c("Average number (across all simulations) of susceptible at end of simulation, as reported by DSAIDE") out_types = rep("Rounded_Integer",nrec) out_notes = rep("Report the rounded integer",nrec) outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task ######################### # Task 4 ######################### tid = tid + 1 tasktext = "Sometimes you do want to know the overall average, including both settings where an outbreak happened and where it didn't. However, often combining those two very different scenarios can be confusing. So let's change our setup to make sure we get outbreaks every time. You learned previously that if there are more initial infected, chances of an outbreak are larger. Let's do that. Keep all settings as before, but change initial number of symptomatic infected to 2. You'll see that there are a few more outbreaks than with 1 infected, but still not all. Keep increasing the starting value for _I_ until all 10 simulations lead to outbreaks." nrec = 1 # number of items to record out_records = c("Minimum number infected such that all 10 simulations produce outbreaks") out_types = rep("Integer",nrec) out_notes = rep("Report the integer",nrec) outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task ######################### # Task 5 ######################### tid = tid + 1 tasktext = "Set initial number infected to 10, set gI=1, leave all other settings as in the previous task (bE=0, bI=0.005, gE=4, w=0.2). Run the simulation for 5 years. You should find an average of around 133 infected at the end. While there is no guarantee that the average of the stochastic model agrees with an equivalent deterministic model for the types of (nonlinear) models we are exploring, they are often close. We can check that by computing the steady states for the deterministic version of the model. If you want to practice doing steady states, you can go ahead and solve the SEIR ODE model for the steady states of the variables. See e.g. the _Patterns of ID_ app for an explanation on how to do it. If you don't want to practice, I'll tell you that the steady state value for S is S=(g~E~ * g~I~)/(b~E~ * g~I~ + b~I~ * g~E~). Use the current values to compute S and compare with results from the simulation." nrec = 2 # number of items to record out_records = c("Theoretical value for susceptible at steady state, based on equation", "Average number (across all simulations) of susceptible at end of simulation") out_types = rep("Rounded_Integer",nrec) out_notes = rep("Report the rounded integer",nrec) outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task ######################### # Task 6 ######################### tid = tid + 1 tasktext = "Let's repeat the previous simulation, but turn off waning immunity. Also set back to gI=2. Still 10 initial infected, random seed at 123, do 20 simulations (be patient), run for 5 years. Record the average number of susceptible left at the end of the simulation. Use that value to compute R~0~ using the final size equation. Then compare it with the R~0~ value computed based on the model (this is the same as for the SIR model, the E compartment does not make a difference here)." nrec = 2 # number of items to record out_records = c("Average number (across all 20 simulations) of susceptible at end of simulation", "R0 based on model equation") out_types = c("Rounded_Integer","Numeric") out_notes = c("Report the rounded integer","Report to one decimal place") outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task ######################### # Task 7 ######################### tid = tid + 1 tasktext = "Finally, let's explore extinctions. We'll start with 100 individuals in each compartment, b~E~=0 and b~I~=0.01, g~E~=4, g~I~=2, w=1. No births and deaths, start simulation at 0, run for 60 months, time step arbitrary. Run 10 simulations with seed of 123. You should see some initial ups and downs, and then for each simulation, you'll see fluctuations around some type of steady state, with an average number of infected at the end of 57. For those settings, in the absence of transmission from the _E_ class, the basic reproductive number (if everyone were susceptible) is R~0~ = S~0~b~I~/g~I~ = 0.01*400/2 = 2. Now, reduce transmission to bI=0.006. Compute R~0~ for this value. Run simulations for 60 months. You will find that for some of these simulations, the pathogen goes extinct (and susceptible go up to the maximum value of 400). Now double simulation time to 10 years. You will find that for most simulations, extinction occured." nrec = 2 # number of items to record out_records = c("Number of simulations that DO NOT lead to extinction after 10 years of simulation (with bI=0.006)", "Reproductive number for bI=0.006") out_types = c("Rounded_Integer","Numeric") out_notes = c("Report the rounded integer","Report to one decimal place") outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task ######################### # Task 8 ######################### tid = tid + 1 tasktext = "Keep everything as in the previous task, but change infected tranmission rate to b~I~=0.007. Compute R~0~ for this value. Run 10 simulations for 10 years. You will find that for most simulations, extinction did not occur." nrec = 2 # number of items to record out_records = c("Number of simulations that DO lead to extinction after 10 years of simulation (with bI=0.007)", "Reproductive number for bI=0.007") out_types = c("Rounded_Integer","Numeric") out_notes = c("Report the rounded integer","Report to one decimal place") outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task ######################### # Task 9 ######################### tid = tid + 1 tasktext = "Keep exploring. You can test how R~0~ and other model quantities affect extinction, and how the longer you run the simulations, the more extinctions you get. You can also turn on births and deaths and see how they affect the overall patterns you find. There is a lot more to explore. For some settings, i.e. if you want to run many replicates for a long time, things are starting to slow down. That's a disadvantage of stochastic models, they often take longer to run." nrec = 1 # number of items to record out_records = c("Nothing") out_types = rep("None",nrec) out_notes = c("") outcomes = data.frame( TaskID = rep(tid,nrec), TaskText = rep(tasktext,nrec), RecordID = paste0('T',tid,'R',(1:nrec)), Record = out_records, Type = out_types, Note = out_notes) alloutcomes = rbind(alloutcomes,outcomes) rc = rc + nrec #increment record counter by number of outcomes to record for this task
#save the fully filled task table to a tsv file alloutcomes$QuizID = paste0(packagename,"_",appsettings$appid) alloutcomes$AppTitle = appsettings$apptitle alloutcomes$AppID = appsettings$appid #remove a few variables from the data frame savedoutcomes <- dplyr::select(alloutcomes,QuizID,AppID,AppTitle,TaskID,TaskText,RecordID,Record,Type,Note) write.table(savedoutcomes, paste0(appsettings$appid,"_tasktable.tsv"), append = FALSE, sep = "\t", row.names = F, col.names = TRUE)
# Take all the text stored in the table and print the tasks and items to record write_tasktext(alloutcomes)
Further Information {#shinytab4}
- This app (and all others) are structured such that the Shiny part (the graphical interface you see and the server-side function that goes with it) calls an underlying R script (or several) which runs the simulation for the model of interest and returns the results.
- For this app, the underlying function(s) running the simulation are called
r appsettings$simfunction. You can call them directly, without going through the shiny app. Use thehelp()command for more information on how to use the functions directly. If you go that route, you need to use the results returned from this function and produce useful output (such as a plot) yourself. - You can also download all simulator functions and modify them for your own purposes. Of course to modify these functions, you'll need to do some coding.
- For examples on using the simulators directly and how to modify them, read the package vignette by typing
vignette('DSAIDE')into the R console. - For some further reading concerning stochastic dynamics of ID, see e.g. [@lloyd-smith05; @black1966; @keeling1997].
References
Try the DSAIDE package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.