---
title: Model Comparison and Norovirus Fit
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 demonstrates basic fitting of data to 2 simple infection models. This shows the concept of model/hypothesis testing. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab. 

__Note:__ You should complete the "Fitting influenza data" app before starting this app. 


## The Model {#shinytab2}

### Model Overview
This app fits 2 different SIR models to norovirus infection data.


#### Models

The overall model is a variant of the basic SIR model, with the inclusion of a process that allows infection of individuals from some common (unmodeled) source.

### Model Diagram
The diagram illustrates the model.

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



### Model Equations
Implementing the models as continuous-time, deterministic systems leads to the following set of ordinary differential equations: 

$$
\begin{aligned}
\dot S & =  - nS  - bSI \\
\dot I &  =  nS + bSI - gI \\
\dot R &  = gI \\
\end{aligned}
$$

### Model Variants
* We explore the following 3 models/hypotheses: 
  * The outbreak involved only person-to-person transmission, i.e. _n=0_.
  * The outbreak involved person-to-person transmission and an outside source of infection (e.g. contaminated food) that lead to constant infections at a fixed rate _n>0_.
  * The outbreak involved person-to-person transmission and an outside source of infection which was present for some period of time during the outbreak, i.e. _n>0_ between some time _t~1~_ and _t~2~_ and 0 otherwise.

### Data source
The data being used in this app is daily new cases of norovirus for an outbreak at a school camp. See `help('norodata') for more details.


### Model comparison 
There are different ways to evaluate how well a model fits to data, and to compare between different models. This app shows the approach of using Akaike's "An Information Criterion" (AIC), or more precisely, we'll use the one corrected for small sample size, AICc . If we fit by minimizing the sum of squares (SSR), as we do here, the formula for the AICc is:
$$
AICc = N \log(\frac{SSR}{N})+2(K+1)+\frac{2(K+1)(K+2)}{N-K}
$$
where _N_ is the number of data points, _SSR_ is the sum of squares residual at the final fit, and _K_ is the number of parameters being fit. A lower value means a better model. One nice feature of the AIC is that one can compare as many models as one wants without having issues with p-values and correcting for multiple comparison, and the models do not need to be nested (i.e. each smaller model being contained inside the larger models). That said, AIC has its drawbacks. Nowadays, if one has enough data available, the best approach is to evaluate model performance by a method like cross-validation [@hastie11].

For evaluation of models with different AIC, there is fortunately no arbitrary, "magic" value (like p=0.05). A rule of thumb is that if models differ by AIC of more than 2, the one with the smaller one is considered statistically better supported (don't use the word 'significant' since that is usually associated with a p-value<0.05, which we don't have here). I tend to be more conservative and want AIC differences to be at least >10 before I'm willing to favor a given model. Also, I think that visual inspection of the fits is useful. If one model has a lower AIC, but the fit doesn't look that convincing biologically (e.g. very steep increases or decreases in some quantity), I'd be careful drawing very strong conclusions.

Note that the absolute value of the AIC is unimportant and varies from dataset to dataset. Only relative differences matter. And it goes without saying that we can only compare models that are fit to exactly the same data.


## 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). If any quantity is not given in those units, you need to convert it first (e.g. if it says a week, you need to convert it to 7 days).*


### Task 1 
* Take a look at the inputs and outputs for the app. It's similar to the "Fitting influenza data" app (which you should do before this one).
* Each model variant fits parameters _b_ and _g_. Model variant 2 also fits a common source infection rate _n_. Model variant 3 additionally fits times _t~1~_ and _t~2~_ at which infection from the common source starts and stops (i.e., is larger than 0). The best fit estimates are shown under the figure, together with the SSR and AICc. 
* For simplicity, lower and upper bounds are set inside the simulator function to the beginning and end of the data. As a result, they can't be adjusted through the user interface.
* To find out more about the data, see _Further Resources_ or `help(norodata)`. 

### Task 2 
* There were a total of 288 susceptible individuals. Use that as starting value for the susceptibles. Assume 1 infected, no recovered.
* Try fit model variant 1, which fits parameters  _b_ and _g_. There are two general approaches. You can set _b_ and _g_ to pretty much any starting value and try to run a large number of iterations with the 3 different solvers until you get a good fit.
* Alternatively, you can try to manually find good starting values for _b_ and _g_ by trying different ones and running a single iteration. Once you have values for which the model is somewhat close to the data, increase iteration, and do the formal fit.
* In theory, you should always end up with the same best fit. However, as discussed in the flu fitting app, optimizers can sometimes get stuck, so it might be for certain settings you don't end up with a good fit.
* The best fit I was able to find for model variant 1 was an SSR=1218.81, AIC = 57.03 and best fit parameter values _b = 0.01_ and _g = 1.21_.  

### Task 3
* Now switch to model variant 2, which also fits parameter _n_.
* Play around with different starting values for the parameters, different optimizers and different numbers of iterations and see what the best fit is you can find.
* The best fit I was able to find is identical to model 1, i.e., the optimizer sets the rate of infection from an external source to _n=0_. This, of course, gives the same SSR. The AIC is now larger at 61.74 since the model has more parameters and gets penalized for this. Based on these 2 models so far, we would conclude that an external source did not improve model fit and only direct transmission (model variant 1) was more likely.

### Task 4
* Now switch to model variant 3, which also fits parameters _t~1~_ and _t~2~_. See what best fit you can find. The value for _t~1~_ must be 8 or greater.    
* I was able, with 5000 iterations of solver 3 and good starting conditions, to get an SSR=17.50 and AIC=34.40  with best fit parameter values _b_/_g_/_n_/_t~1~_/_t~2~_ of 0.0079/1.6357/0.9960/10.7444/12.0000. You might be able to find an even better fit (lower SSR).

The SSR is clearly lower than that for model 1 and 2, and the AIC value suggests that the extra model parameters are 'worth it', and overall, this model should be favored.

You might have noticed that getting a good fit is tricky and often you don't reach a good one. This comes back to the 'optimizer getting stuck' concept discussed in the 'flu fitting' 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.
* The data for this study is saved in the data variable `norodata`, you can read more about it by looking at its help file entry `help(norodata)`. The publication from which the data comes is [@kuo09].

### References


