Exploring the impact of parameter changes

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 allows you to explore the effect of specific model parameters on some outcomes of interest for the simple SIR model with births and deaths that we explored previously. 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.

Note that the model we use here is the same as the one in the Stochastic SIRS model app, but you don't have to have worked through that app yet to follow along, you will have seen all the model components in other apps previously.

Learning Objectives

• Learn how to use models to explore the impact of parameters on outcomes

The Model {#shinytab2}

Model Overview

The model used here is the SIR model with births and deaths and waning immunity. This model tracks susceptibles, infected/infectious and recovered hosts. The following compartments are included:

• S - uninfected and susceptible individuals
• I - individuals who are infected and infectious.
• R - recovered/removed individuals. Those individuals have recovered and are immune.

The included processes/mechanisms are the following:

• Susceptible individuals can become infected at rate b.
• Infected hosts recover at rate g.
• Immunity can wane at rate w.
• New susceptible hosts enter the system (are born) at rate n. Natural mortality 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 model implemented as set of differential equations is given by the following equations:

$$\dot S = n - bSI - mS + wR$$ $$\dot I = bSI - gI - mI$$ $$\dot R = gI - mR - wR$$

Model exploration

The new component implemented for this app is the ability to run the simulation multiple times as you scan through values for a specific parameter. Instead of running the simulation once for a given choice of parameter values and looking at the resulting time-series, this app lets you explore the impact of each model parameter on some pre-specified outcomes. Specifically, we consider the maximum and final value of each variable.

In the app, those outcomes of interest are labeled Smax, Imax, Rmax and Sfinal, Ifinal, Rfinal. You can choose one of the model parameters to be varied between some minimum and maximum value. The other parameter values remain fixed. For each parameter value, the model is run and the outcomes of interest computed. The resulting plot is one showing how those outcomes of interest vary with the parameter you investigated.

While you can do that with some of the other apps too, by manually changing parameters, re-running the simulation, and recording results, this is done here in an automated way. Exploring models in this way, and presenting plots showing how some outcome(s) of interest vary with a specific parameter, is a very common way models are used in research.

What to do {#shinytab3}

The tasks below are described in a way that assumes everything is in units of days (rate parameters, therefore, have units of inverse days).

#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 initial number of susceptible/infected/recovered to 1000/1/0. Set infection rate to 0.002, and recovery rate to 1. The former value does not matter since we'll be exploring the simulation for different values of _b_. Set natural births and deaths to 0 and no waning immunity. Choose _b_ as the parameter to vary, go from 0.0005 to 0.005, do it for 10 different parameter values, linear spacing, and no log scales for plotting. Set start and final simulation time to 0 and 100, time steps 0.1. Run the simulation. You will notice that you get a message saying the system did not reach steady state 1 times (see output message below plot). Increase simulation time to 1000 days, now you should reach steady state for all 10 simulations. Recall that for a simple SIR model with no births/deaths or waning immunity, the final steady state is one in which the outbreak is over and possibly some individuals remain uninfected. 

In this app, you don't see the model dynamics, i.e. the time-series for the variables. It is being run in the background, but only parts of it, i.e. maximum and final states of variables are recorded. Try to visualize what the S/I/R curves look like, based on your exploration of the simple SIR model app. Then check these 'virtual' curves with the results reported in the plot and make sure you understand how each of them comes about. For instance do you understand why _Smax_ is always 1000? Or why _Rfinal_ increases from 0 to 1000 as you increase _b_? If you have a hard time with this, I suggest you (re)visit the basic SIR model and explore it some more. The final number of susceptibles declines as transmission (and thus the reproductive number) increases. This hopefully all makes sense to you based on what you have learned in the other apps. The maximum and final values for I and R should equally make sense. Essentially, you are running a simulation 10 times that can produce (at most) a single outbreak. If transmission is low enough, no outbreak occurs. The larger the transmission rate, the larger the outbreak. Instead of showing the time-series for each outbreak, the simulation only records some summary quantities and reports those. This is a very common way of analyzing models and reporting findings. Switch the plot engine to plotly and use it to read off the value for the maximum number infected for b=0.002, you should find it to be 154. The read off the maxium number infected for b=0.005.

Play around with the sample number, i.e. number of different parameter values (the higher the number, the more often the simulation runs and the longer it takes), the minimum and maximum, and linear or logarithmic spacing and plotting, to understand how these different settings impact the resulting plot."

nrec = 1 # number of items to record
out_records = c("Maximum number infected for b = 0.005")
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 = "A good way to build intuition about a specific model is to run simulations, look at the results, and then look at the model equations and diagram and see if by staring at them, you can understand why you got the results you saw. With enough practice, it is often possible to intuit specific results based on the underlying equations - but that requires a good bit of modeling practice. Let's try to comparison of model equations and simulation. Look at the basic model with no births/deaths/waning immunity and try to predict how changing the parameter _g_ will impact the reported outcomes. You can also think of it scientifically, namely as _g_ increases, the duration of the infectious period becomes shorter, what does that mean for the outbreak? Now explore this with the model. Set the parameter to change to be _g_, values from 0.2 to 2. Fix b=0.001. Run 10 simulations for 1000 days. Try both linear and logarithmic spacing. Take a look at the results and see if they agree with what you expected. Use the plotly plot again to report the requested value." 

nrec = 1 # number of items to record
out_records = c("Final number susceptible for g=1")
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 = "In the previous task, you should have noticed that above _g=1_ the final number of recovered (which is also the max) goes to 0. What does that mean for the outbreak? From the previous apps, you know that the criterion for getting an outbreak or not is R0=1. Let's compute R0 here for this model (with no births/deaths/waning immunity). If you don't remember how to do this, re-visit the apps in the _Reproductive Number_ section. Then also compute the R0 for the smallest and largest values for _g_ in the ranges you explored in the previous task."

nrec = 1 # number of items to record
out_records = c("Reproductive number for g=0.2")
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, the model is simple enough that one can figure things out just doing math, without needing to run simulations. For the simple model we have here, we can compute the steady states for the different variables. This was one of the tasks in the _Patterns of ID_ app. Revisit that app if you don't remember. We found there that at steady state, in the presence of births and deaths (and ignoring waning immunity), the variable S is defined by model parameters through  _S = (g+m)/b_, and similar equations for the other variables. Based on this, we expect that the final, steady state value for S increases linearly with _g_. Let's compare this equation with running the simulation. Set n=10, m = 0.01 (these values will keep the population at a steady value of 1000 in the absence if infection. If you are unclear why that is so, revisit the _ID Patterns_ app). Keep everything as in the previous task, i.e. we are still scanning over g from 0.2 to 2 with b=0.001. Compare the values for _S~final~_ for the different values of _g_ you get from the mathematical equation with what you get from running the simulation. Things should be consistent between the math and the simulations. If they are not, it means something went wrong somewhere. You should see that S~final~ increases linearly with _g_, up to the point where R0 drops below one, thus no outbreak occurs and the steady state in the presence of infection (which the equation we use here specifies) is not valid anymore. You can test the other parts of the equation by varying _m_ (the equations tells you that you should again see a linear change) and _b_ (you should see an inverse relationship between increasing _b_ and _S_).

Even for this simple model, we cannot easily compute mathematical expressions for the maximum values of the variables. In general, as soon as our model reaches a certain level of complexity (maybe around 5 equations and more), getting exact analytic/mathematical equations for most outcomes of interest is not possible. Also, as models get more complicated, even experienced modelers can often not intuit what model behavior one should expect as specific parameters vary. Often the only way to find out is by actually running the simulations. "

nrec = 1 # number of items to record
out_records = c("Final number INFECTED for g=0.4")
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 

#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)

As an end note, it is worth pointing out that for this app, we focus on varying a single parameter at a time and looking at the results. With some coding effort, we could change the underlying simulation to loop over say 2 parameters and produce outcomes for sets of parameter values. The results could be plotted as a 2-dimensional heatmap for each outcome. While this could be extended to more than 2 parameters, it will become hard to visualize and will also take longer to run. What is often done is that some parameters are chosen at discrete values (e.g. g=0.1 and 1) and then for those sets, one loops over some other parameter (e.g. b). This can be done in any combinations that are deemed interesting to explore. If there are many parameters that one might want to change, a different approach is useful, which you'll learn about in the Uncertainty and Sensitivity app.

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 running the simulation is 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.

References

DSAIDE documentation built on Aug. 24, 2023, 1:07 a.m.