---
title: "Simulating disease progress curves"
author: "Kaique S Alves"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simulating disease progress curves}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## Introduction

In `epiffiter`, opposite to the `fit_` functions (estimate parameters from fitting models to the data), the `sim_` family of functions allows to produce the DPC data given a set of parameters for a specific model. Currently, the same four population dynamic models that are fitted to the data can be simulated.

The functions use the `ode()` function of the `devolve` package (Soetaert,Petzoldt & Setzer 2010) to solve the differential equation form of the e epidemiological models.

## Hands On

### Packages

First, we need to load the packages we'll need for this tutorial.

```{r message=FALSE, warning=FALSE}
library(epifitter)
library(magrittr)
library(ggplot2)
library(cowplot)
```

### The basics of the simulation

The `sim_` functions, regardless of the model, require the same set of six arguments. By default, at least two arguments are required (the others have default values)

-   `r`: apparent infection rate
-   `n`: number of replicates

When `n` is greater than one, replicated epidemics (e.g. replicated treatments) are produced and a level of noise (experimental error) should be set in the `alpha` argument. These two arguments combined set will generate `random_y` values, which will vary randomly across the defined number of replicates.

The other arguments are:

-   `N`: epidemic duration in time units
-   `dt`: time (fixed) in units between two assessments
-   `y0`: initial inoculum
-   `alpha`: noise parameters for the replicates

#### Exponential

Let's simulate a curve resembling the exponential growth.

```{r}
exp_model <- sim_exponential(
  N = 100,    # total time units 
  y0 = 0.01,  # initial inoculum
  dt = 10,    #  interval between assessments in time units
  r = 0.045,  #  apparent infection rate
  alpha = 0.2,# level of noise
  n = 7       # number of replicates
)
head(exp_model)
```

A `data.frame` object is produced with four columns:

-   `replicates`: the curve with the respective ID number
-   `time`: the assessment time
-   `y`: the simulated proportion of disease intensity
-   `random_y`: randomly simulated proportion disease intensity based on the noise

Use the [`ggplot2`](https://ggplot2.tidyverse.org/) package to build impressive graphics!

```{r}
exp_plot = exp_model %>%
  ggplot(aes(time, y)) +
  geom_jitter(aes(time, random_y), size = 3,color = "gray", width = .1) +
  geom_line(size = 1) +
  theme_minimal_hgrid() +
  ylim(0,1)+
  labs(
    title = "Exponential",
    y = "Disease intensity",
    x = "Time"
  )
exp_plot
```

#### Monomolecular

The logic is exactly the same here.

```{r}
mono_model <- sim_monomolecular(
  N = 100,
  y0 = 0.01,
  dt = 5,
  r = 0.05,
  alpha = 0.2,
  n = 7
)
head(mono_model)
```

```{r}
mono_plot = mono_model %>%
  ggplot(aes(time, y)) +
  geom_jitter(aes(time, random_y), size = 3, color = "gray", width = .1) +
  geom_line(size = 1) +
  theme_minimal_hgrid() +
  labs(
    title = "Monomolecular",
    y = "Disease intensity",
    x = "Time"
  )
mono_plot
```

#### The Logistic model

```{r}
logist_model <- sim_logistic(
  N = 100,
  y0 = 0.01,
  dt = 5,
  r = 0.1,
  alpha = 0.2,
  n = 7
)
head(logist_model)
```

```{r}
logist_plot = logist_model %>%
  ggplot(aes(time, y)) +
  geom_jitter(aes(time, random_y), size = 3,color = "gray", width = .1) +
  geom_line(size = 1) +
  theme_minimal_hgrid() +
  labs(
    title = "Logistic",
    y = "Disease intensity",
    x = "Time"
  )
logist_plot
```

#### Gompertz

```{r}
gomp_model <- sim_gompertz(
  N = 100,
  y0 = 0.01,
  dt = 5,
  r = 0.07,
  alpha = 0.2,
  n = 7
)
head(gomp_model)
```

```{r}
gomp_plot = gomp_model %>%
  ggplot(aes(time, y)) +
  geom_jitter(aes(time, random_y), size = 3,color = "gray", width = .1) +
  geom_line(size = 1) +
  theme_minimal_hgrid() +
  labs(
    title = "Gompertz",
    y = "Disease intensity",
    x = "Time"
  )
gomp_plot
```

### Combo

Use the function `plot_grid()` from the [`cowplot`](https://wilkelab.org/cowplot/index.html) package to gather all plots into a grid

```{r fig.height=6, fig.width=8}
plot_grid(exp_plot,
          mono_plot,
          logist_plot,
          gomp_plot)
```

## References

Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer (2010). Solving Differential Equations in R: Package deSolve. Journal of Statistical Software, 33(9), 1--25. DOI: [10.18637/jss.v033.i09](http://dx.doi.org/10.18637/jss.v033.i09)