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

[pwinskill]

Our previous implementation of a flexible carrying capacity that the user could set using set_ carrying_capacity() asked for absolute values of the carrying capacity at time change points to be specified. This in turn required the user to know the baseline carrying capacity (k0), which required set_equilibrium() to be called. This was problematic for downstream processes such as calibration. When calibrating we want to be able to re-set the equilibrium (by re-calling set_equilibrium()) whilst retaining the impact of any previously set interventions. As set_ carrying_capacity() previously depended on k0 this broke the work flow.

To address this the user now specifies the carrying capacity changes as scalers which rescale k0 and are therefore indepent of the aboslute value of k0. This change retains all previous options:

# Illustration that new method retains the flexibility of the old

# Baseline run
timesteps <- 365 * 2
p1 <- malariasimulation::get_parameters(
  overrides = list(
    human_population = 1000,
    progress_bar = TRUE
  )
) |>
  malariasimulation::set_equilibrium(init_EIR = 10)

set.seed(1)
s1 <- malariasimulation::run_simulation(
  timesteps = timesteps,
  parameters = p1
)

# Decreasing carrying capacity at timepoint
p2 <- p1 |>
  set_carrying_capacity(timesteps = 180, carrying_capacity_scalers = matrix(0.01))

set.seed(1)
s2 <- malariasimulation::run_simulation(
  timesteps = timesteps,
  parameters = p2
)

# Increasing  carrying capacity at timepoint
p3 <- p1 |>
  set_carrying_capacity(timesteps = 180, carrying_capacity_scalers = matrix(5))

set.seed(1)
s3 <- malariasimulation::run_simulation(
  timesteps = timesteps,
  parameters = p3
)

# Decreasing then increasing then original  carrying capacity at timepoints
p4 <- p1 |>
  set_carrying_capacity(timesteps = c(180, 180 + 150, 180 + 300), carrying_capacity_scalers = matrix(c(0.01, 5, 1)))

set.seed(1)
s4 <- malariasimulation::run_simulation(
  timesteps = timesteps,
  parameters = p4
)

# Illustration that the impact is now maintained after re-calling set_equilibrium()
# Impact of setting carrying capacity is maintained when re-setting init EIR
p5 <- p2 |>
  set_equilibrium(init_EIR = 1)

set.seed(1)
s5 <- malariasimulation::run_simulation(
  timesteps = timesteps,
  parameters = p5
)

plot(s1$timestep, s1$n_detect_730_3650 / s1$n_730_3650, t = "l", ylim = c(0, 1), ylab = "PfPr")
lines(s2$timestep, s2$n_detect_730_3650 / s2$n_730_3650, col = "green")
lines(s3$timestep, s3$n_detect_730_3650 / s3$n_730_3650, col = "red")
lines(s4$timestep, s4$n_detect_730_3650 / s4$n_730_3650, col = "blue")
lines(s5$timestep, s5$n_detect_730_3650 / s5$n_730_3650, col = "orange")
abline(v = 180, lty = 2)

The new method also is species specific:

# Can still get different scalers for different species
timesteps <- 365 * 2
p1 <- malariasimulation::get_parameters(
  overrides = list(
    human_population = 1000,
    progress_bar = TRUE
  )
) |>
  malariasimulation::set_equilibrium(init_EIR = 10)

p6 <- p1 |>
  set_species(
    species = list(arab_params, fun_params, gamb_params),
    proportions = c(0.335, 0.335, 0.33)
  ) |>
  malariasimulation::set_equilibrium(init_EIR = 10) |>
  set_carrying_capacity(timesteps = 180, carrying_capacity_scalers = matrix(c(0.1, 1, 2), ncol = 3))

set.seed(1)
s6 <- malariasimulation::run_simulation(
  timesteps = timesteps,
  parameters = p6
)

plot(s6$total_M_arab ~ s6$timestep, t = "l", ylim = c(0, 5000), ylab = "total_M")
lines(s6$total_M_fun ~ s6$timestep, col = "red")
lines(s6$total_M_gamb ~ s6$timestep, col = "blue")

The new method also can be combined with seasonal profiles:

# Can still set seasonality and carrying capacity rescaling
p7 <- malariasimulation::get_parameters(
  overrides = list(
    human_population = 1000,
    progress_bar = TRUE,
    model_seasonality = TRUE,
    g0 = 0.284596,
    g = c(-0.317878,-0.0017527,0.116455),
    h = c(-0.331361,0.293128,-0.0617547)
  )
) |>
  set_carrying_capacity(timesteps = c(365, 365 * 2), carrying_capacity_scalers = matrix(c(0.5, 0.25))) |>
  malariasimulation::set_equilibrium(init_EIR = 10)

set.seed(1)
s7 <- malariasimulation::run_simulation(
  timesteps = timesteps * 2,
  parameters = p7
)

plot(s7$timestep, s7$n_detect_730_3650 / s7$n_730_3650, t = "l", ylim = c(0, 1), ylab = "PfPr")  

[RJSheppard]

All looks good to me! I wonder whether it would be good to reference and explain how the initial carrying capacity gets generated at the start of the vignette. i.e., you don't explicitly set carrying capacity parameter, but the carrying capacity gets generated in set_equilibrium, to match the EIR and population size, and this is what the scalers are relative to. (Assuming I've understood this correctly!)

[tbreweric]

Really nice! I had one suggestion about wording in the markdown but otherwise looks great!

[giovannic]

If we're sure we'll never use absolute carrying capacity here, then this is cleaner and feel free to merge.

If not, then it would probably be best to copy set_carrying_capacity to set_relative_carrying_capacity which sets the appropriate parameters.

Thanks for this!

[giovannic]

Also, to pass the merge checks please run devtools::document() to update the documentation and push up the generated files

[pwinskill]

Thanks all! 🙏😎

@giovannic - I've left get_init_carrying_capacity() as a user-facing function. So in the case where someone would want to work with absolute carrying capacity they could just convert in to relative before specifiying. I can't imagine that it will be very common.