---
title: Model Comparison Fit
output:
  html_document:
    theme: null
    highlight: null
    fig_retina: null
    fig_caption: true
    mathjax: default 
    keep_md: false
bibliography: references.bib
---

```{r, include = FALSE}
#*************************************
#general setup to define package and get path locations
#all paths are inside the package and retrieved with system.file
packagename = "DSAIRM"
helperdir = "helperfunctions"
mbmodeldir = "mbmodels"
figuredir = "media"
appdocdir = "appinformation" 
#*************************************
#Note: for this to process/knit, several helper functions need to be available (sourced) first
#those are in the inst/helperfunctions folder
#Note: in general, the "processing-script.R" in the docsfordevelopers should be used to produce the html docs
#manual knitting of each doc only during development/testing
#*************************************
files_to_source = list.files(system.file(helperdir,package = packagename),full.names=TRUE)
sapply(files_to_source, source) #sourcing needs to happen inside each Rmd file since knitr starts a new environment
#load the settings file for the current app 
#so we can automatically include figure, list the functions in the further information section
#and use other information specific to the current app for the task table generation
currentrmdfile = knitr::current_input() 
appsettings = get_settings(currentrmdfile,appdocdir,packagename)
```


## Overview {#shinytab1}
This app shows some 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.


## The Model {#shinytab2}

### Model Overview
This app fits 2 different models to virus load data from humans.

#### Model 1

This model consists of 4 compartments. The following entities are modeled:

* **U** - uninfected cells 
* **I** - infected cells
* **V** - (free) virus
* **X** - T-cells

The following processes/flows are included in the model: 

* Virus infects cells at rate _b_.
* Infected cells produce new virus at rate _p_, are killed by T-cells at rate _k_ and die due to other causes at rate _d~I~_.
* Free virus is removed at rate _d~V~_ or goes on to infect further uninfected cells. It can also infect new cells at rate _b_, with unit conversion factor _g_.
* T-cells are activated proportional to virus at rate _a_ and undergo exponential growth at rate _r_.


#### Model 2

This model also consists of 4 compartments. The following entities are modeled:

* **U** - uninfected cells 
* **I** - infected cells
* **V** - (free) virus
* **X** - B-cells/antibodies

The following processes/flows are included in the model: 

* Virus infects cells at rate _b_.
* Infected cells produce new virus at rate _p_ and die at rate _d~I~_.
* Free virus is removed by antibodies at rate _k_, and by other mechanisms at rate _d~V~_. It can also infect new cells at rate _b_, with unit conversion factor _g_.
* B-cells/antibodies grow exponentially proportional to virus at rate _a_ and decay at rate _d~X~_.


### Model Diagram
The diagram illustrates both compartmental models, with colors to indicate mechanisms that are part of either model 1 or model 2.

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



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

$$\dot U =  - bUV$$ 
$$\dot I = bUV - d_I I - kIX$$
$$\dot V = pI - d_V V - gbUV$$
$$\dot X = a V + r X$$


Model 2 is described by these equations:
$$\dot U =  - bUV$$ 
$$\dot I = bUV - d_I I $$
$$\dot V = pI - d_V V - k V X - gbUV$$
$$\dot X = a V X - d_X X$$

### Comments
Note that both models are rather simple and - arguably - somewhat unrealistic. For instance the T-cells grow exponentially without any limit. This model can be justified if one only tries to fit data from an acute infection - as done here for influenza - where the infection is over (i.e. virus drops below level of detection) before T-cells reach their peak. For both models, we ignored birth/death of uninfected cells, which is possibly ok for an acute infection (but needs to be determined based on known biology for a given system). 

We also assume for model 2 that the dynamics of B-cells and antibodies track each other close enough that we can lump them together into a single equation.

In general, models need to be carefully tailored to the question and setting. This is true for models used in fitting even more than for models that are used without trying to perform data fitting.


### Data source
The data being used in this app comes from [@hayden96]. Specifically, the data being fit is the 'no intervention' viral load data shown in Figure 2 of the paper.

### Fitting Approach
The basic fitting approach is the one described in the _Basic model fitting_ app, i.e. we minimize the sum of squares of the log differences between model predictions and virus load data, taking censored values into account.


### 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 "Basic Fitting" app (which you should do before this one).
* Each model has 3 parameters that are being estimated, the best fit estimates are shown under the figure, together with the SSR and AICc. Switch to model 2, run and compare how it differs from model 1.

### Task 2:
* You probably noticed that the default setting is only 10 iterations for the solver. That's way too few to find the best fit, so comparing the 2 models after such few iterations is not very meaningful. Increase the iterations to 100 or more and re-run both models.
* You'll likely find that even with 100 iterations, the fits don't look that great. Explore different starting values for the parameters that are being fit (a,b and r/dX) and see if you can get better fits. You can also explore different values for the fixed parameters. Usually, those should not be changed since they are fixed based in biological a priori information. If we didn't know those parameter values, we would also try to fit them - though that can lead to overfitting, as you learned in another of the fitting apps. For this exercise, you can play with those fixed values.

### Task 3:
* Based on the AICc for both models after 100 or more iterations, what do you conclude regarding the quality of each model? 
* Take a look at SSR for both models. You should see the model with the lower AICc also have the lower SSR. Why is that so? Is that always true for any 2 models? (It's not: The 2 models here are in some sense special, understand how/why.)
* If you haven't figured out the last question, and also for a general check, compute AICc 'by hand' for each of the models using the reported SSR and the formula given in the "Models" section. You should of course find the values shown below the plot.
* By doing all this, you should have realized that since the models have the same numbers of parameters (and are fit to the same data, which always has to be the case if you want to compare models), all the terms in the AICc equation other than SSR are the same for the 2 models. Thus, in this case, SSR and AICc behave the same. If one of the models had more parameters than the other, you could find a model with lower SSR but higher AIC (because the additional complexity of the model, evidenced by more parameters, doesn't lower the SSR enough for the AIC to go down.)


## 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('DSAIRM')` into the R console.

### References


