Merge pull request #272 from mrc-ide/carrying_capacity_relative

Custom carrying capacity is set as a relative scaler, not absolute value

R/compartmental.R

@@ -16,7 +16,7 @@ parameterise_mosquito_models <- function(parameters, timesteps) {
         for(j in 1:length(parameters$carrying_capacity_timesteps)){
           timeseries_push(
             k_timeseries,
-            parameters$carrying_capacity_values[j,i],
+            parameters$carrying_capacity_scalers[j,i] * k0,
             parameters$carrying_capacity_timesteps[j]
           )
         }

R/variables.R

@@ -1,4 +1,5 @@
-#' @title Define model variables #' @description
+#' @title Define model variables 
+#' @description
 #' create_variables creates the human and mosquito variables for
 #' the model. Variables are used to track real data for each individual over
 #' time, they are read and updated by processes

R/vector_control_parameters.R

@@ -143,23 +143,24 @@ set_spraying <- function(
 #' 
 #' @param parameters the model parameters
 #' @param timesteps vector of timesteps for each rescale change
-#' @param carrying_capacity matrix of baseline carrying_capacity for each species 
-#' With nrows = length(timesteps), ncols = length(species)
+#' @param carrying_capacity_scalers matrix of scaling factors to scale the baseline 
+#' carrying capacity for each species with nrows = length(timesteps),
+#'  ncols = length(species)
 #' 
 #' @export
 set_carrying_capacity <- function(
     parameters,
     timesteps,
-    carrying_capacity
+    carrying_capacity_scalers
 ){
-  stopifnot(nrow(carrying_capacity) == length(timesteps))
-  stopifnot(ncol(carrying_capacity) == length(parameters$species))
+  stopifnot(nrow(carrying_capacity_scalers) == length(timesteps))
+  stopifnot(ncol(carrying_capacity_scalers) == length(parameters$species))
   stopifnot(min(timesteps) > 0)
-  stopifnot(min(carrying_capacity) >= 0)
+  stopifnot(min(carrying_capacity_scalers) >= 0)
 
   parameters$carrying_capacity <- TRUE
   parameters$carrying_capacity_timesteps <- timesteps
-  parameters$carrying_capacity_values <- carrying_capacity
+  parameters$carrying_capacity_scalers <- carrying_capacity_scalers
   parameters
 }
 

src/RcppExports.cpp

@@ -260,7 +260,7 @@ BEGIN_RCPP
 END_RCPP
 }
 
-RcppExport SEXP run_testthat_tests();
+RcppExport SEXP run_testthat_tests(void);
 
 static const R_CallMethodDef CallEntries[] = {
     {"_malariasimulation_create_adult_mosquito_model", (DL_FUNC) &_malariasimulation_create_adult_mosquito_model, 5},

tests/testthat/test-vector-control.R

@@ -372,18 +372,18 @@ test_that('set_carrying_capacity works',{
       species = "gamb",
       carrying_capacity = TRUE,
       carrying_capacity_timesteps = 1,
-      carrying_capacity_values = matrix(0.1)
+      carrying_capacity_scalers = matrix(0.1)
     )
   )
 
   expect_error(
     set_carrying_capacity(p, 1, matrix(c(0.1, 0.1), nrow = 2)),
-    "nrow(carrying_capacity) == length(timesteps) is not TRUE",
+    "nrow(carrying_capacity_scalers) == length(timesteps) is not TRUE",
     fixed = TRUE
   )
   expect_error(
     set_carrying_capacity(p, 1, matrix(c(0.1, 0.1), ncol = 2)),
-    "ncol(carrying_capacity) == length(parameters$species) is not TRUE",
+    "ncol(carrying_capacity_scalers) == length(parameters$species) is not TRUE",
     fixed = TRUE
   )
   expect_error(
@@ -393,7 +393,7 @@ test_that('set_carrying_capacity works',{
   )
   expect_error(
     set_carrying_capacity(p, 1, matrix(-1)),
-    "min(carrying_capacity) >= 0",
+    "min(carrying_capacity_scalers) >= 0",
     fixed = TRUE
   )
 })

vignettes/Carrying-capacity.Rmd

@@ -81,17 +81,13 @@ lines(s_seasonal$pfpr ~ s_seasonal$timestep, col = "darkorchid3")
 ## Custom modification of the mosquito carrying capacity
 
 In malariasimulation, we can also modify the carrying capacity over time in a more
-bespoke manner. Firstly it is helpful to start with our baseline (equilibrium)
-carrying capacity. We can access this for our baseline parameters with:
+bespoke manner. We can use the `set_carrying_capacity()` helper function to 
+specify scalars that modify the species-specific carrying capacity through time,
+relative to the baseline carrying capacity, at defined time steps. The baseline
+carrying capacity is generated when calling `set_equilibrium()`, to match the 
+specified EIR and population size.
 
-```{r}
-cc <- get_init_carrying_capacity(p)
-cc
-```
-
-We can then use this in combination with `set_carrying_capacity()` to provide
-modified, species-specific carrying capacity values at given timepoints. This
-allows us to model a number of useful things:
+This functionality allows us to model a number of useful things:
 
 ### Larval source management
 
@@ -102,10 +98,10 @@ can be used to specify this intervention
 ```{r, fig.width = 7, fig.height = 4, out.width='100%'}
 # Specify the LSM coverage
 lsm_coverage <- 0.8
-# Set LSM by reducing the carrying capacity by (1 - coverage)
+# Set LSM by scaling the carrying capacity by (1 - coverage)
 p_lsm <- p |>
   set_carrying_capacity(
-    carrying_capacity = matrix(cc * (1 - lsm_coverage), ncol = 1),
+    carrying_capacity_scalers = matrix(1 - lsm_coverage, ncol = 1),
     timesteps = 365
   )
 set.seed(123)
@@ -124,7 +120,7 @@ abline(v = 365, lty = 2, col = "grey60")
 ### Invasive species
 
 An invasive species may exploit a new niche, increase the carrying
-capacity at the point of invasion. The functions
+capacity at the point of invasion. The function
 `set_carrying_capacity()` gives complete freedom to allow
 changes to the carrying capacity to be made. Here
 we will demonstrate using carrying capacity rescaling to capture the
@@ -136,13 +132,13 @@ invasion of _Anopheles stephensi_.
 p_stephensi <- p |>
   set_species(list(gamb_params, steph_params), c(0.995, 1 - 0.995)) |>
   set_equilibrium(init_EIR = 5)
-cc_invasive <- get_init_carrying_capacity(p_stephensi)
+
 # Next, at the time point of invasion, we scale up the carrying capacity for
-# the invasive species by a factor that will be dependent on the proporties of
+# the invasive species by a factor that will be dependent on the properties of
 # the invasion.
 p_stephensi <- p_stephensi |>
   set_carrying_capacity(
-    carrying_capacity = matrix(cc_invasive * c(1, 2000), ncol = 2),
+    carrying_capacity_scalers = matrix(c(1, 2000), ncol = 2),
     timesteps = 365
   )
 
@@ -169,7 +165,7 @@ time
 ```{r, fig.width = 7, fig.height = 4, out.width='100%'}
 p_flexible <- p |>
   set_carrying_capacity(
-    carrying_capacity = cc * matrix(c(0.1, 2, 0.1, 0.5, 0.9), ncol = 1),
+    carrying_capacity = matrix(c(0.1, 2, 0.1, 0.5, 0.9), ncol = 1),
     timesteps = seq(100, 500, 100)
   )