#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)
This app allows you to explore a simple SIR model with 2 types of hosts. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.
This app assumes that you have worked through several of the previous ones, e.g. that you are familiar with the reproductive number and how it is computed.
This model tracks susceptibles, infected and recovered of 2 different types. Think of those types as e.g. males/females, children/adults, etc.
The following compartments are included, twice for each type (i=1,2):
The included processes/mechanisms are the following:
The flow diagram and equations describe the model implemented in this app:
knitr::include_graphics(here::here('inst/media',appsettings$modelfigname))
$$\dot S_1 = - S_1 (b_{11} I_1 + b_{12} I_2) + w_1 R_1 $$ $$\dot I_1 = S_1 (b_{11} I_1 + b_{12} I_2) - g_1 I_1 $$ $$\dot R_1 = g_1 I_1 - w_1 R_1 $$ $$\dot S_2 = - S_2 (b_{21} I_1 + b_{22} I_2) + w_2 R_2 $$ $$\dot I_2 = S_2 (b_{21} I_1 + b_{22} I_2) - g_2 I_2 $$ $$\dot R_2 = g_2 I_2 - w_2 R_2 $$
It might be worth saying something about the transmission terms. I generally use a single subscript to describe transmission from a group, e.g. b~A~ for transmission/infectiousness of asymptomatic. If there are multiple groups that can be susceptible and infectious, a common notation is to start with the receiving group first, then the sending/transmitting group, e.g. if susceptible individuals of type 1 are infected by individuals of type 2, most authors write b~12~. I follow this convention. Note however that it is equally ok to use b~12~ to mean that infected type 1 individuals transmit to susceptible type 2 individuals. I actually like this sender first perspective/notation better, and used it originally, but switched to stick with the convention used in the main introductory textbooks on this topic.
In general you need to read papers/model descriptions carefully, and hopefully the authors do a good job explaining exactly what is meant. Such that there is no confusion. Just read carefully every time and don't jump to conclusions based on what you have seen before or what you think it means.
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).
#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 = " Start with 1000 susceptible hosts and 1 infected host of type one. 200 susceptible hosts and 1 infected host of type two. Simulation duration approximately 5 years. Assume that transmission from host 1 to host 1 is _b~11~ = 0.002_, from host 2 to host 2 is _b~22~ = 0.01_. No transmission from one host type to the other _b~12~ = 0_ and _b~21~ = 0_. Assume that the duration of the infectious period is 1 month long for both types of hosts (i.e. same recovery rate). No waning immunity. Run the simulation and ensure you get outbreaks in both populations with approximately 20% susceptibles left at the end." nrec = 2 # number of items to record out_records = c("Number susceptible of type 1 left at end of simulation", "Number susceptible of type 2 left 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 = "Set _b~22~ = 0.02_. Rest as before. Run the simulation. You should get the same outbreak as before among type 1 hosts, a larger outbreak among type 2 hosts. For our current choice of parameters, more specifically the transmission rates, it is ok to define and compute separate R~0~ for the two populations. Contemplate why it is ok to do so and compute the two R~0~ values. Use what you learned about R~0~ to compute a theoretical value, then check with the simulation using the final size equation." nrec = 2 # number of items to record out_records = c("R0 for type 1 hosts", "R0 for type 2 hosts") 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 3 ######################### tid = tid + 1 tasktext = "Set _b~11~ = 0.001_. Rest as in previous task. Run the simulation. You should get the same outbreak as before among type 2 hosts, but no real outbreak among type 1 hosts (though a few infections will occur). Compute the R~0~ values for these settings and convince yourself that the theoretical values and simulation results agree." nrec = 2 # number of items to record out_records = c("R0 for type 1 hosts", "R0 for type 2 hosts") 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 = "Now set the transmission rate to host 1 from host 2 _b~12~ = 0.001_. Everything else unchanged. Run the simulation. You should see an outbreak in both populations. This is an example of a (small) core group driving an outbreak in the larger group. Now that the two groups interact, the simple individual equations for _R~0~_ do not apply anymore. One can compute an overall _R~0~_ for the joint populations, but that is a bit involved and beyond what we want to do here. For more details on that, see e.g. [@keeling07]." nrec = 2 # number of items to record out_records = c("Number susceptible of type 1 left at end of simulation", "Number susceptible of type 2 left 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 5 ######################### tid = tid + 1 tasktext = "Let's test the idea that the small core group is driving the outbreaks and it's not just because we introduced cross-transmission and that feature alone led to outbreaks in both populations. To do so, set _b~12~_ to 0 and instead set _b~21~ = 0.001_. If it was just due to cross-transmission, we should probably see outbreaks in both populations. If it was a more infectious/transmissible (i.e. higher _R~0~_) core-group driving the outbreak in the less transmissible group, we should not see an outbreak. Run the simulation and interpret the results." 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 ######################### # Task 6 ######################### tid = tid + 1 tasktext = "Now set transmission rates to _b~11~_ = _b~22~_ = 0 and _b~12~_ = 0.01, _b~21~_ = 0.002. Note that these are the same values for the transmission terms as in task 1, but we now there is no transmission among individuals of the same type and non-zero transmission between types. Contemplate what you expect to see, run the simulation, see if your expectations are confirmed. Compare the results to those from task 1." nrec = 2 # number of items to record out_records = c("Number susceptible of type 1 left at end of simulation", "Number susceptible of type 2 left 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 7 ######################### tid = tid + 1 tasktext = "Note that things are not symmetric here when it comes to values of cross-transmission. Explore this by switching the values for the cross-transmission terms, such that now _b~12~_ = 0.002, _b~21~_ = 0.01. Leave everything else as in the previous task. Run the simulation and compare the results to those from the previous task." nrec = 2 # number of items to record out_records = c("Number susceptible of type 1 left at end of simulation", "Number susceptible of type 2 left 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 8 ######################### tid = tid + 1 tasktext = "You have previously encountered the idea that transmission occurs only between hosts of different types, namely in the vector-borne transmission app. There we thought of the different types as different species, e.g. humans and mosquitoes. But it could also be among different types of the same hosts. For instance this type of transmission could represent a sexually transmitted disease in a heterosexual population, with the 2 types of hosts being females and males. Many sexually transmitted infections do not produce life-long immunity, re-infection is possible and not uncommon. Let's explore this. Set the number of susceptibles for both populations to 1000, 1 infected in each population. Set b~12~ = 0.002, b~21~ = 0.001. Set the other transmission terms to 0. Leave recovery rates as before. Then turn on waning immunity, assume it has an average duration of 5 months for each population (i.e. rates _w~1~_ and _w~2~_ need to be the inverse of 5 months). Run the simulation for 10 years. Confirm that you reach a steady state, with around 667 susceptible of type 1 at the end." nrec = 1 # number of items to record out_records = c("Number susceptible of type 2 left 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 9 ######################### tid = tid + 1 tasktext = "In the previous task, the steady state values for S/I/R for the 2 populations were different. Think about why that is the case. Try to find the value for b~21~ for which the 2 populations reach the same steady state. Then explore different values for both transmission rates and see how that impact the outcomes. Consider what that means for a real infectious disease. Would you expect transmission rates to be the same or not? What does it depend on?" nrec = 1 # number of items to record out_records = c("Value for b21 at which steady state values for both populations are the same.") out_types = rep("Numeric",nrec) out_notes = rep("Report the exact value",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 10 ######################### tid = tid + 1 tasktext = "Keep exploring. You can continue the exploration from the previous task by altering other parameters for the 2 types. You might for instance want to consider some real ID where accounting for 2 types of hosts is important, e.g. some type of STI. You can try to go to the literature and see if you can find parameter values for the ID of interest and explore potential patterns of that ID with the simulation." 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
### #These tasks are not really doing what I expect them to do, not sure what's going on just now, no time to figure out ### # tasktext = "You learned that the reproductive number is defined as new infections/infectious hosts of the same type produced by one infectious host. For the case where this is cross-transmission, one can derive it without too much effort. This follows the same reasoning you saw in the vector-borne transmission app. Type 1 infected produce new infections of type 2 at rate _b~21~_ for time 1/g~1~, and if initially all type 2 are susceptible, this leads to a total number of new type 2 infected of _b~21~ S~2~_/_g~1~_. A similar equations holds for the number of type 1 infections produced by a type 2 infected host. The total reproductive number is the product, thus _R~0~_ = _S~1~ S~2~_ _b~21~_ _b~12~_ / (_g~1~ g~2~_). Compute _R~0~_ for the current model settings." # tasktext = "Let's explore the idea of something like a sexually transmitted disease that only goes from host 1 to 2 and the other way (assuming a fully heterosexual population) and combine it with the core group idea. Above, we found that the reproductive number was the product of the 2 components. Let's compute them separately. You should find that for the settings above, they are equal and both larger than 1. Let's change that. Set the 2 populations to equal size, both 1000 susceptible. Set _b~12~_ = 0.01, _b~21~_ = 0.002. Now compute both the expected new infections of type 2 produced by type 1, R~{21}, and the equivalent value for R~{12}, as well as the overall _R~0~_. What do you find? Given those values, do you expect to see an outbreak when you run the simulation with those values? Run and confirm."
#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)
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. vignette('DSAIDE')
into the R console.Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.