---
title: "koalas"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{koalas}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=8, fig.height=4
)

options("pboptions" = list(type = "none"))
```

## Overview

The koalas package is intended to be used for exploring chlamydia control programmes.  The package is still being developed, and feedback is welcomed (contact details are available [here](https://ivh.ku.dk/english/employees/?pure=en/persons/487288) and via my private consultancy company [here](http://www.m2d2consultancy.com)).

## Installation

The latest release version of the package can be installed from r-universe:

```{r eval=FALSE}
install.packages(c("koalas","tidyverse"),
  repos=c("https://ku-awdc.r-universe.dev/", "https://cran.rstudio.com/"))
```

This also installs the tidyverse package from CRAN, which is used in this vignette for plotting etc. Note that C++ compilers are NOT required for installation as long as your installed version of R is current. This should cover most use cases.

Alternatively, you can install the development version directly from GitHub.  To do that you must first pre-install the remotes package (from CRAN) and the IPDMR package (from ku-awdc.r-universe.dev/):

```{r eval=FALSE}
install.packages(c("IPDMR","remotes"),
  repos=c("https://ku-awdc.r-universe.dev/", "https://cran.rstudio.com/"))
```

Then you can install the koalas package from source:

```{r eval=FALSE}
remotes::install_github("ku-awdc/koalas")
```

Note that this last step requires C++ compilers.


## TLDR; Quick guide to running scenarios

The pre-determined scenarios can be run in the following way:

```{r}
library("koalas")
model <- KoalasV2$new()
model$burnin()
```

This runs the model from 1st July 2022 up to 1st October 2025, which is where we might start doing active interventions:

```{r}
library("ggplot2")
theme_set(theme_light())

model$autoplot()
```

### Baseline scenario

For the baseline scenario, we can just clone the model and then run it for a further 10 years with no interventions:

```{r}
baseline <- model$clone(deep=TRUE)
baseline$run(10, frequency=0)
baseline$autoplot() +
  geom_vline(xintercept=baseline$run_dates, lty="dashed")
```

### Testing interventions (phase 1)

For other scenarios, we split into 2 phases:

1. Phase 1: high-intensity sampling
2. Phase 2: lower-intensity sampling

For each of these we can use the assess_interventions function - for example:

```{r}
phase1 <- assess_interventions(
  model,
  years = 2,
  frequency = 1:4,
  prop_active = seq(0, 1, by=0.1),  # Change to by=0.01 for more precise results
  prop_targeted = 0.0
)

library("dplyr")
phase1 |>
  filter(Prevalence <= 5) |>
  group_by(Frequency) |>
  arrange(PropActive) |>
  slice(1L)
```

NOTE: these are rough numbers - change to by=0.01 for more precise results!

Then we can plot these results:

```{r}
phase1 |>
  ggplot(aes(x=PropActive, y=Prevalence, col=Frequency)) +
  geom_line() +
  geom_hline(yintercept=5, lty="dashed") +
  theme(legend.position="bottom")
```

We can also look at the final koala populations:

```{r}
phase1 |>
  ggplot(aes(x=PropActive, y=Koalas, col=Frequency)) +
  geom_line() +
  theme(legend.position="bottom")
```

Or we can look at the effect of targeted interventions alone:

```{r}
phase1_targeted <- assess_interventions(
  model,
  years = 2,
  frequency = 1:4,
  prop_active = 0.0,
  prop_targeted = seq(0, 1, by=0.1)  # Change to by=0.01 for more precise results
)
phase1_targeted |>
  ggplot(aes(x=PropTargeted, y=Prevalence, col=Frequency)) +
  geom_line() +
  geom_hline(yintercept=5, lty="dashed") +
  theme(legend.position="bottom") +
  ylim(0,100)
```

Or we can look at a range of active proportions on top of targeted proportions fixed at 80%:

```{r}
phase1_both <- assess_interventions(
  model,
  years = 2,
  frequency = 1:4,
  prop_active = seq(0, 1, by=0.1),  # Change to by=0.01 for more precise results
  prop_targeted = 0.8
)
phase1_both |>
  ggplot(aes(x=PropActive, y=Prevalence, col=Frequency)) +
  geom_line() +
  geom_hline(yintercept=5, lty="dashed") +
  theme(legend.position="bottom") +
  ylim(0,100)
```


### Testing interventions (phase 2)

For phase 2, we first need to update the model using a phase 1 strategy e.g. a frequency of 4 and sampling proportion of 47% (the minimum required to get prevalence under 5% after 2 years):

```{r}
model$run(2, frequency = 4, prop_active = 0.47, prop_targeted = 0.0)
```

We can see that this does reduce the prevalence to less than 5%:

```{r}
model$autoplot()
```

Then we can assess phase 2 interventions using this as a starting point:

```{r}
phase2 <- assess_interventions(
  model,
  years = 8,
  frequency = 1:4,
  prop_active = seq(0, 1, by=0.1),  # Change to by=0.01 for more precise results
  prop_targeted = 0.0
)
phase2 |>
  filter(Prevalence <= 5) |>
  group_by(Frequency) |>
  arrange(PropActive) |>
  slice(1L)
```

NOTE: these are rough numbers - change to by=0.01 for more precise results!

Again, we can plot these results (a different graph is shown here):

```{r}
library("tidyr")

phase2 |>
  pivot_longer("Prevalence":"Koalas", names_to="Metric", values_to="Value") |>
  ggplot(aes(x=PropActive, y=Value, col=Frequency)) +
  geom_line() +
  facet_wrap(~Metric, ncol=1, scales="free_y") +
  geom_hline(data=tibble(Metric="Prevalence", ll=5), mapping=aes(yintercept=ll), lty="dashed")
```

### Final scenario

Finally, we can update the model for 8 years with a frequency of 2 and proportion of 45% (the minimum required to keep the prevalence under 5%):

```{r}
model$run(8, frequency = 2, prop_active = 0.45)
model$autoplot() +
  geom_vline(xintercept=model$run_dates, lty="dashed")
```


## Creating a model object

The koalas model is designed as an encapsulated object, which is created as follows:

```{r}
library("koalas")
model <- KoalasV2$new()
```

This creates a model that you can then interrogate and manipulate using the active bindings and methods supplied. For example, we can set the initial state:

```{r}
model$set_state(S=299, I=1, R=0, Af=0, Cf=0)
```

And one or more parameters:

```{r}
model$set_parameters(birthrate=0.38)
```

For a complete list of valid state and parameter values see the help file for the class:  ?KoalasV2

You can also see the current state and parameters using their active bindings:

```{r}
model$state |> simplify2array()
model$parameters |> simplify2array()
```


## Updating the model

To update the model for one or more day, use the update method with given number of days:

```{r}
model$update(21)
```

You can then extract the results so far as data frame:

```{r}
model$results_wide
```

When you want to apply an active intervention, use the active_intervention method with specified proportion of animals to test/treat/vaccinate:

```{r}
model$active_intervention(proportion=0.5)
```

You can then run for another period of time and re-extract results:

```{r}
model$update(21)
model$results_wide
```

Note that results are cumulative, i.e. you can keep updating the model and it will carry on from where it left off. This means you can apply interventions and/or change parameter values at whatever time points are required.

To re-set the model, create a new instance:

```{r}
model <- KoalasV2$new()
model$update(5)
model$results_wide
```

You can also extract aggregated results in a format more suitable for plotting:

```{r}
model$results_long
```

This produces the following (non-exclusive) combinations:

- Healthy = S+V+N+R
- Infectious = I+If+Af+Cf
- Diseased = Af+Cf
- Infertile = Sf+Vf+Nf+Rf+I+Af+Cf
- Immune = V+Vf+R+Rf+N+Nf




## Package versions

```{r}
sessionInfo()
```