Running_EPIONCHO_IBM.Rmd

---
title: "Installing and Running EPIONCHO-IBM"
author:
- name: Jonathan Hamley
  affiliation: Imperial College London, University of Bern
  role:         # Contributorship roles (e.g., CRediT)
  - Conceptualization
  - Writing - Original Draft Preparation
  - Writing - Review & Editing
  - Software
- name: Martin Walker
  affiliation: Royal Veterinary College, Imperial College London
  role:         
  - Conceptualization
  - Software
- name: Philip Milton
  affiliation: Imperial College London
  role:         
  - Conceptualization
  - Writing - Original Draft Preparation
  - Writing - Review & Editing
  - Software
- name: Matthew A Dixon
  affiliation: Imperial College London, SCI Foundation 
  role:         
  - Software
  - Data/ code curation 
date: "Febuary 14, 2023"
output:
  html_document:
    df_print: paged
  word_document: default
  radix::radix_article: default
  pdf_document: default
description: |
  A stochastic, individual-based model for onchocerciasis transmission
---

**EPIONCH0-IBM is a stochastic, individual-based model for onchocerciasis transmission. A mathematical description can be found at https://doi.org/10.1371/journal.pntd.0007557.**

# 1. System requirements

EPIONCHO IBM is written in R (1) which is freely available for installation on Windows, Mac OS and Linux platforms from https://www.r-project.org. It is also recommended to install R Studio https://www.rstudio.com which provides a helpful user interface for R. EPIONCHO-IBM has been tested on Mac OS and Windows, but should be compatible with Linux and run on any desktop or laptop machine.

# 2. Installation guide

The model package can be downloaded and installed directly from a GitHub repository (https://github.com/mrc-ide/EPIONCHO.IBM). The remotes (2) package must be installed to do this. Installation should take no more than a few seconds on most desktop or laptop machines.


```{r, eval=T, echo = FALSE, collapse= FALSE}
remotes::install_github("mrc-ide/EPIONCHO.IBM")
library(EPIONCHO.IBM)
```

# 2.1 Running the model for multiple runs

To see how to deal with multiple runs, and processing data over multiple runs, the Ov16 vignette has an example of running the model multiple times, along with a function to process the data from multiple runs.

# 3. Demo

## 3.1 Simulating an endemic equilibrium without MDA

Before simulating treatment, the model must be simulated to the endemic equilibrium. It is best to run the model to equilibrium (```run_equilibrium = TRUE```, ```give.treat = 0 ```), save the output and use this as the starting point for MDA. Note that running the model to the equilibrium can be slow, taking in excess of 5 minutes (depending on machine performance). To run the model to equilibrium without treating in the same run, use the following code:

```{r}
#length of simulation in years
timesteps = 30

#should treatment be given, when and how often
give.treat.in = 0
treat.strt = 1; treat.stp = 16
trt.int = 1

#annual biting rate, which determines infection prevalence (60% microfilarae prevalence)
ABR.in = 1082

output_equilibrium <- ep.equi.sim(time.its = timesteps,
                                  ABR = ABR.in,
                                  N.in = 440,
                                  treat.int = trt.int,
                                  treat.prob = 0.65,
                                  give.treat = give.treat.in,
                                  treat.start = treat.strt,
                                  treat.stop = treat.stp,
                                  pnc = 0.05,
                                  min.mont.age = 5,
                                  vector.control.strt = NA,
                                  gam.dis.in = 0.3,
                                  run_equilibrium = TRUE,
                                  equilibrium = NA,
                                  print_progress = TRUE)
```

To stop the printing of the time as the model runs, set ```print_progress = FALSE```. Currently, with ```print_progress = TRUE```, when the model run reaches each year in the simulation, the year and % progress will be printed.
If ```run_equilibrium = TRUE``` then no input is expected for ```equilibrium```, which is where a previously simulated endemic equilibrium can be entered. This will be discussed later in section 3.3. Above, no MDA is not requested by setting ```gv.trt = 0```. Note that values must be entered for treatment start, stop and interval (1 is annual, 0.5 is biannual), e.g ```treat.strt = 0; treat.stp = 16; trt.int = 1```, even if ```gv.trt = 0```. 

It is also possible to manually change the density-dependent parameters relating to the establishment of the adult *Onchocerca volvulus* (```delta.hz.in```, ```delta.hinf.in``` and ```c.h.in```) and individual exposure heterogeneity in humans (```gam.dis.in```), in the function calling the model (above). However, if you use a individual exposure heterogeneity of 0.2, 0.3, or 0.4, the density dependent parameters will be set automatically, to their respectively fitted density-dependent parameters in [Hamley et al. 2019](https://journals.plos.org/plosntds/article?id=10.1371/journal.pntd.0007557). 

## 3.2 Visualising microfilarial infection prevalence through time

The output from run_EPIONCHO_IBM is a list with 21 elements.

```{r}
names(output_equilibrium)
```

Within the first 10 elements, there are 5 main outputs (`mf_prev`, `mf_intens`, `ov16_seroprevalence_no_seroreversion`, `ov16_seroprevalence_finite_seroreversion`, `L3`, and `blackfly_l3_prevalence`) which contain the temporal dynamics (at 1 day increments) for the microfilarial prevalence (individuals age > 5), the population mean microfilarial intensity (individuals age > 5), the temporal dynamics for seroprevalence (w/o and w/ seroreverstion), the mean number of L3 larvae in the blackfly population, and the estimated prevalence of l3 in the blackfly population. To plot these over time, use: 

```{r}
years <- output_equilibrium$years

plot(years, output_equilibrium$mf_prev, type = 'l', xlab = 'time (years)', ylab = 'microfilarial prevalence', ylim = c(0, 1))

plot(years, output_equilibrium$mf_intens, type = 'l', xlab = 'time (years)', ylab = 'mean microfilarial intensity')

plot(years, output_equilibrium$ov16_seroprevalence_no_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ no (black) and w/ finite (red) seroreversion', ylim = c(0, 1))
lines(years, output_equilibrium$ov16_seroprevalence_finite_seroreversion, col="red")

plot(years, output_equilibrium$L3, type = 'l', xlab = 'time (years)', ylab = 'mean L3 in blackfly population')

plot(years, output_equilibrium$blackfly_l3_prevalence, type = 'l', xlab = 'time (years)', ylab = 'L3 prevalence in blackfly population')

# worm burden parameter explained below
plot(years, output_equilibrium$worm_burden_outputs[, "fertile_female_worm_burden"], type = 'l', xlab = 'time (years)', ylab = 'fertile female worm burden')
```

The `ABR` attribute lets you know what the input ABR used was. 
The `all_infection_burdens` attribute is a matrix that contains the detailed infection burden of each individual at the final timestep.
The `ov16_indiv_matrix` attribute is a matrix/dataframe that contains the detailed ov16 serostatus for each individual at various timesteps. See the Ov16 vignette for more detail.
The `years` attribute lets you know what the timestep is in terms of years from timepoint 0.
The `ABR_recorded` attribute is a temporal output showing the ABR at every year/timestep.
The `coverage.recorded` attribute is a temporal output showing the MDA coverage at every year/timestep. It will have NA values if no treatment was given.
The `percent_never_treated` attribute is a temporal output, showing how many individuals at a given timestep had never been treated. Without MDA, this will be all NA values.
The last element,```"all_equilibrium_outputs"``` is the element required to simulate MDA without having to run the model to the endemic equilibrium.

The `all_mf_prevalence_age_grouped`, `all_mf_intensity_age_grouped`, `ov16_timetrend_outputs`, `ov16_timetrend_outputs_adj`, `worm_burden_outputs`, `oae_incidence_outputs`, and `all_morbidity_prevalence_outputs` elements contain the temporal dynamics (at 1 day increments) for multiple age groups for mf prevalence, intensity, ov16 seroprevalence (w/ and without seroreversion), ov16 seroprevalence adjusted, worm burden (male, fertile female, infertile female), OAE incidence, and morbidity prevalence (including the OAE prevalence) respectively.
The prevalence and burden outputs are grouped by [5, 80], [0-5), [0-5), [5-10), [10-15), [15-20), [20-30), [30-50), [50-80].
Other than 5-80, the rest of the age groups can be customized by editing the `output_age_groups` input parameter.
The OAE incidence outputs are grouped by  [5, 80], [0-5), [5-10], [11-15], Males, and Females. These cannot be customized. 

## 3.3 Simulating annual and biannual mass drug administration with ivermectin
### 3.3.1 Annual MDA

Assuming you have already run the model and saved the output as described in section 3.1, you can simulate 15 years of annual mass drug administration with ivermectin by running:

```{r}
timesteps = 30 
treat.strt = 1; treat.stp = 26 #if model run stops in year 26, the last round is at the beginning of year 25
gv.trt = 1
trt.int = 1

output_treat_annual <- ep.equi.sim(time.its = timesteps,
                                   ABR = ABR.in,
                                   N.in = 440,
                                   treat.timing = NA,
                                   treat.int = trt.int,
                                   treat.prob = 0.65,
                                   give.treat = gv.trt,
                                   treat.start = treat.strt,
                                   treat.stop = treat.stp,
                                   pnc = 0.05,
                                   min.mont.age = 5,
                                   vector.control.strt = NA,
                                   gam.dis.in = 0.3,
                                   run_equilibrium = FALSE,
                                   equilibrium = output_equilibrium$all_equilibrium_outputs,
                                   print_progress = FALSE)

```

If in the above model run ```print_progress = FALSE``` is called, but ```gv.trt = 1``` is also specified, only the duration of MDA will be printed (no progress updates). Currently, with ```print_progress = TRUE```, when the model run reaches each year in the simulation, the year and % progress will be printed as before.
Note that ```timesteps``` is the total number of years for which the model is run. ```treat.start``` and ```treat.stop``` are the years within this duration at which MDA is given.

To visualise the infection dynamics through time during annual MDA, use:

```{r}
years <- output_treat_annual$years

plot(years, output_treat_annual$mf_prev, type = 'l', xlab = 'time', ylab = 'microfilarial prevalence', ylim = c(0, 1))

abline(v = seq(1, 25), col = 'grey', lwd = 0.1)

plot(years, output_treat_annual$mf_intens, type = 'l', xlab = 'time', ylab = 'mean microfilarial intensity')

abline(v = seq(1, 25), col = 'grey', lwd = 0.1)

plot(years, output_treat_annual$ov16_seroprevalence_no_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ no (black) and w/ finite (red) seroreversion', ylim = c(0, 1))
lines(years, output_treat_annual$ov16_seroprevalence_finite_seroreversion, col="red")

abline(v = seq(1, 25), col = 'grey', lwd = 0.1)

plot(years, output_treat_annual$percent_never_treated, type = 'l', xlab = 'time', ylab = 'Percent never treated')

abline(v = seq(1, 25), col = 'grey', lwd = 0.1)

plot(years, output_treat_annual$worm_burden_outputs[, "fertile_female_worm_burden"], type = 'l', xlab = 'time (years)', ylab = 'fertile female worm burden')
abline(v = seq(1, 25), col = 'grey', lwd = 0.1)
```

### 3.3.2 Biannual MDA

Assuming you have run the model to the endemic equilibrium (as described in secton 3.1), temporal infection trends during biannual MDA can be obtained by running:

```{r}
timesteps = 30 
treat.strt = 1; treat.stp = 26 #if treatment stops in year 26, the last round is at the beginning of year 25
gv.trt = 1
trt.int = 0.5


output_treat_biannual <- ep.equi.sim(time.its = timesteps,
                                   ABR = ABR.in,
                                   N.in = 440,
                                   treat.timing = NA,
                                   treat.int = trt.int,
                                   treat.prob = 0.65,
                                   give.treat = gv.trt,
                                   treat.start = treat.strt,
                                   treat.stop = treat.stp,
                                   pnc = 0.05,
                                   min.mont.age = 5,
                                   vector.control.strt = NA,
                                   gam.dis.in = 0.3,
                                   run_equilibrium = FALSE,
                                   equilibrium = output_equilibrium$all_equilibrium_outputs,
                                   print_progress = TRUE)
```

Note than above ```trt.int = 0.5``` (rather than ```trt.int = 1```) to request biannual MDA. The output can be visualised using:

```{r}
years <- output_treat_biannual$years

plot(years, output_treat_biannual$mf_prev, type = 'l', xlab = 'time', ylab = 'microfilarial prevalence', ylim = c(0, 1))

abline(v = seq(1, 25, 0.5), col = 'grey', lwd = 0.1)

plot(years, output_treat_biannual$mf_intens, type = 'l', xlab = 'time', ylab = 'mean microfilarial intensity')

abline(v = seq(1, 25, 0.5), col = 'grey', lwd = 0.1)

plot(years, output_treat_biannual$ov16_seroprevalence_no_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ no (black) and w/ finite (red) seroreversion', ylim = c(0, 1))
lines(years, output_treat_biannual$ov16_seroprevalence_finite_seroreversion, col="red")

plot(years, output_treat_biannual$percent_never_treated, type = 'l', xlab = 'time', ylab = 'Percent never treated')

abline(v = seq(1, 25), col = 'grey', lwd = 0.1)

plot(years, output_treat_biannual$worm_burden_outputs[, "fertile_female_worm_burden"], type = 'l', xlab = 'time (years)', ylab = 'fertile female worm burden')
lines(years, output_treat_biannual$worm_burden_outputs[, "male_worm_burden"], type = "l", col = "red")
lines(years, output_treat_biannual$worm_burden_outputs[, "infertile_female_worm_burden"], type = "l", col = "blue")
abline(v = seq(1, 25), col = 'grey', lwd = 0.1)
```

### 3.3.3 Coverage and non-compliance

The population level coverage is a proportion and is entered via ```treat.prob```. Note that because children <5 do not receive treatment, coverage of 1 cannot be achieved. Additionally, there is the option to set the proportion of individuals who never take treatment (```pnc = 0.05```), or using a correlation parameter for a more dynamic compliance structure (```comp.correlation = 0.3```). More information on correlation can be seen in the complex MDA vignette.

### 4. References

1. R Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, 2018).

2. Csárdi, G. RCurl: Package ‘remotes’: R Package Installation from Remote Repositories, Including 'GitHub'. (2022). https://cran.r-project.org/web/packages/remotes/