Installation
The latest release version of the package can be installed from
r-universe:
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/):
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:
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:
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:
library("ggplot2")
#> Warning: package 'ggplot2' was built under R version 4.3.3
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:
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:
- Phase 1: high-intensity sampling
- Phase 2: lower-intensity sampling
For each of these we can use the assess_interventions function - for
example:
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")
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
phase1 |>
filter(Prevalence <= 5) |>
group_by(Frequency) |>
arrange(PropActive) |>
slice(1L)
#> # A tibble: 3 × 5
#> # Groups: Frequency [3]
#> Frequency PropActive PropTargeted Prevalence Koalas
#> <fct> <dbl> <dbl> <dbl> <dbl>
#> 1 2 0.8 0 3.66 182.
#> 2 3 0.6 0 4.26 177.
#> 3 4 0.5 0 3.77 174.
NOTE: these are rough numbers - change to by=0.01 for more precise
results!
Then we can plot these results:
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:
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%:
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):
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%:
Then we can assess phase 2 interventions using this as a starting
point:
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)
#> # A tibble: 4 × 5
#> # Groups: Frequency [4]
#> Frequency PropActive PropTargeted Prevalence Koalas
#> <fct> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0.8 0 1.71 284.
#> 2 2 0.5 0 1.22 256.
#> 3 3 0.4 0 0.290 243.
#> 4 4 0.3 0 0.520 242.
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):
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%):
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:
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:
model$set_state(S=299, I=1, R=0, Af=0, Cf=0)
And one or more parameters:
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:
model$state |> simplify2array()
#> S V I N R Af Cf Sf Vf If Nf Rf Day
#> 299 0 1 0 0 0 0 0 0 0 0 0 0
#> SumTx SumVx SumRx SumMx
#> 0 0 0 0
model$parameters |> simplify2array()
#> vacc_immune_duration vacc_redshed_duration
#> 1.000 0.500
#> natural_immune_duration beta
#> 1.000 2.750
#> subcinical_duration subclinical_recover_proportion
#> 0.500 0.350
#> diseased_recover_proportion birthrate
#> 0.000 0.380
#> acute_duration lifespan_natural
#> 0.400 6.000
#> lifespan_acute lifespan_chronic
#> 4.000 4.800
#> relative_fecundity sensitivity
#> 0.000 0.950
#> specificity cure_prob_N
#> 0.999 0.900
#> cure_prob_I cure_prob_A
#> 0.900 0.750
#> cure_prob_C vaccine_efficacy
#> 0.600 0.500
#> vaccine_booster passive_intervention_rate
#> 0.500 0.030
Updating the model
To update the model for one or more day, use the update method with
given number of days:
You can then extract the results so far as data frame:
model$results_wide
#> # A tibble: 22 × 19
#> Date Day S V I N R Af Cf Sf
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2022-01-01 0 299 0 1 0 0 0 0 0
#> 2 2022-01-02 1 299. 0.0123 1.01 7.47e-6 6.74e-5 4.75e-7 2.29e-14 0
#> 3 2022-01-03 2 299. 0.0246 1.01 1.50e-5 1.37e-4 3.76e-6 1.77e-12 1.79e-18
#> 4 2022-01-04 3 299. 0.0369 1.02 2.25e-5 2.09e-4 1.25e-5 2.13e-11 2.94e-17
#> 5 2022-01-05 4 300. 0.0491 1.02 3.01e-5 2.86e-4 2.93e-5 1.21e-10 2.20e-16
#> 6 2022-01-06 5 300. 0.0614 1.03 3.77e-5 3.70e-4 5.66e-5 4.64e-10 1.05e-15
#> 7 2022-01-07 6 300. 0.0736 1.04 4.53e-5 4.60e-4 9.67e-5 1.38e- 9 3.75e-15
#> 8 2022-01-08 7 300. 0.0858 1.04 5.30e-5 5.59e-4 1.52e-4 3.46e- 9 1.10e-14
#> 9 2022-01-09 8 300. 0.0981 1.05 6.07e-5 6.67e-4 2.24e-4 7.64e- 9 2.79e-14
#> 10 2022-01-10 9 300. 0.110 1.06 6.85e-5 7.86e-4 3.15e-4 1.53e- 8 6.33e-14
#> # ℹ 12 more rows
#> # ℹ 9 more variables: Vf <dbl>, If <dbl>, Nf <dbl>, Rf <dbl>, Z <dbl>,
#> # SumTx <dbl>, SumVx <dbl>, SumRx <dbl>, SumMx <dbl>
When you want to apply an active intervention, use the
active_intervention method with specified proportion of animals to
test/treat/vaccinate:
model$active_intervention(proportion=0.5)
You can then run for another period of time and re-extract
results:
model$update(21)
model$results_wide
#> # A tibble: 43 × 19
#> Date Day S V I N R Af Cf Sf
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2022-01-01 0 299 0 1 0 0 0 0 0
#> 2 2022-01-02 1 299. 0.0123 1.01 7.47e-6 6.74e-5 4.75e-7 2.29e-14 0
#> 3 2022-01-03 2 299. 0.0246 1.01 1.50e-5 1.37e-4 3.76e-6 1.77e-12 1.79e-18
#> 4 2022-01-04 3 299. 0.0369 1.02 2.25e-5 2.09e-4 1.25e-5 2.13e-11 2.94e-17
#> 5 2022-01-05 4 300. 0.0491 1.02 3.01e-5 2.86e-4 2.93e-5 1.21e-10 2.20e-16
#> 6 2022-01-06 5 300. 0.0614 1.03 3.77e-5 3.70e-4 5.66e-5 4.64e-10 1.05e-15
#> 7 2022-01-07 6 300. 0.0736 1.04 4.53e-5 4.60e-4 9.67e-5 1.38e- 9 3.75e-15
#> 8 2022-01-08 7 300. 0.0858 1.04 5.30e-5 5.59e-4 1.52e-4 3.46e- 9 1.10e-14
#> 9 2022-01-09 8 300. 0.0981 1.05 6.07e-5 6.67e-4 2.24e-4 7.64e- 9 2.79e-14
#> 10 2022-01-10 9 300. 0.110 1.06 6.85e-5 7.86e-4 3.15e-4 1.53e- 8 6.33e-14
#> # ℹ 33 more rows
#> # ℹ 9 more variables: Vf <dbl>, If <dbl>, Nf <dbl>, Rf <dbl>, Z <dbl>,
#> # SumTx <dbl>, SumVx <dbl>, SumRx <dbl>, SumMx <dbl>
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:
model <- KoalasV2$new()
model$update(5)
model$results_wide
#> # A tibble: 6 × 19
#> Date Day S V I N R Af Cf Sf Vf If
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2022-01-01 0 300 0 0 0 0 0 0 0 0 0
#> 2 2022-01-02 1 300. 0.0123 0 0 0 0 0 0 0 0
#> 3 2022-01-03 2 300. 0.0247 0 0 0 0 0 0 0 0
#> 4 2022-01-04 3 300. 0.0370 0 0 0 0 0 0 0 0
#> 5 2022-01-05 4 301. 0.0493 0 0 0 0 0 0 0 0
#> 6 2022-01-06 5 301. 0.0616 0 0 0 0 0 0 0 0
#> # ℹ 7 more variables: Nf <dbl>, Rf <dbl>, Z <dbl>, SumTx <dbl>, SumVx <dbl>,
#> # SumRx <dbl>, SumMx <dbl>
You can also extract aggregated results in a format more suitable for
plotting:
model$results_long
#> # A tibble: 36 × 4
#> Date Compartment Koalas Percent
#> <date> <fct> <dbl> <dbl>
#> 1 2022-01-01 Total 300 100
#> 2 2022-01-01 Healthy 300 100
#> 3 2022-01-01 Infectious 0 0
#> 4 2022-01-01 Diseased 0 0
#> 5 2022-01-01 Infertile 0 0
#> 6 2022-01-01 Immune 0 0
#> 7 2022-01-02 Total 300. 100
#> 8 2022-01-02 Healthy 300. 100
#> 9 2022-01-02 Infectious 0 0
#> 10 2022-01-02 Diseased 0 0
#> # ℹ 26 more rows
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