---
title: Exploring the impact of parameter changes 
output:
  html_document:
    theme: null
    highlight: null
    fig_retina: null
    fig_caption: true
    mathjax: default 
    keep_md: false
bibliography: references.bib
---

```{r, echo = FALSE}
#this code loads the settings file for the current app so we can automatically 
#list the functions in the further information section
currentrmdfile = knitr::current_input() 
currentappinfo = gsub("_documentation.Rmd" ,"_settings.R",currentrmdfile)
source(currentappinfo)
```

## Overview {#shinytab1}
This app allows you to explore the effect of specific model parameters on some outcomes of interest for the simple bacteria model. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.


## The Model {#shinytab2}

### Model Overview
The model used here is the SIR model with births and deaths. It is also used and described in the stochastic SIR app. 
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 (S) can become infected at rate _b_. 
* Infected hosts recover at rate _g_. 
* New susceptible hosts enter the system (are born) at rate _m_. Natural death occurs at rate _n_.



### Model Implementation
The flow diagram for the model implemented in this app is:

```{r BIdiagram,  fig.cap='Flow diagram for this model. ',  echo=FALSE}
knitr::include_graphics("../media/stochasticSIRmodel.png")
```


The deterministic model implemented as set of differential equations is given by the following equations:

$$\dot S = m - bSI - nS$$
$$\dot I = bSI - gI - nI$$
$$\dot R = gI - nR$$

This is almost the same model as the basic SIR model from the introductory app, with the only difference that this model also allows natural births and deaths. 


### Model exploration
The new component implemented for this app is the ability to systematically run the simulation multiple times as you scan of 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 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).*

### Task 1 
* Set 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.
* Choose _b_ as the parameter to vary, go from 0.0005 to 0.005, do 10 different parameter values, linear spacing, and no log scales for plotting.
* Set maximum simulation time to 100. Run simulation. You will notice that you get a message saying the system did not reach steady state 10 times, i.e. for each simulation (see output message below plot). Increase simulation time to 1000 days, now you should reach steady state for 5 of the runs.

### Task 2 
* Play around with the 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.

* 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.


### Task 3 
* Another way to build intuition about a given system is by taking a look at the model equations and diagram and see if by 'staring' at them, you can understand why you got the results you saw in the last task. 

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. 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.

### Task 4
The model we have here is so simple that we can in fact figure out some of the outcomes shown in the plot without having to run simulations but instead by doing some math. Specifically, we can compute the steady states for the different variables. This was one of the tasks in the _Patterns of ID_ app. We repeat it here.

Once the system has settled down, there are no more changes in the numbers for each compartment. Mathematically, that means that the left side of the differential equations becomes 0, and they turn into the following algebraic equations: $0 = m - bSI - nS$, $0 = bSI - gI - nI$ and $0 = gI - nR$. One can solve those equations for each of the compartments to get a mathematical expression of what _S_/_I_/_R_ are at steady state. Try to do this. You should find that _S~steady~ = (g+n)/b_ and a similar equations for the other variables. 

### Task 5 
* Compare the mathematical equations you found in the previous task with what you find in the plots. Things should be consistent, e.g. since the mathematical equation says that _S~steady~_ should increase linearly with _g_, you should find that when you run the simulation and look at the plot.

We cannot compute similar 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 analytic/mathematical equations for most outcomes of interest is not possible and the numerical approach of running the simulations and looking at the results is the only option we have.

### Task 6 
* Let's explore the impact of a different parameter. Set all variables as in task 1.
* Choose _g_ as the parameter to vary and set lower and upper limit to 0.4 and 4.
* Before you run the model, try to predict what you should see in the plot for the different outcomes.
* Now run the model, compare the results with your expectations. How do they agree or disagree?

### Task 7 
* In the previous task, you should have noticed that above _g=2_ the final number of recovered (which is also the max) goes to 0. Try to understand what happens (qualitatively) with your model around that value and how that leads to the results you are seeing. (Hint: It has to do with the reproductive number.) 

### Task 8
* Do the same explorations you did above for any other parameter you want to investigate. 

Note that here we focus on a single parameter at a time. 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 long to run. If there are many parameters that could 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


