Parameterise the model using a baseline prevalence (pfpr2-10)

[giovannic]

The current model does this by:

• running a simulation
• executing this function, (with side-effects on the simulation)
• and finding the root of the function outputs here

We would like to create behaviour like:

simparams <- calibrate_model_for_prevalence(simparams, pfpr2_10)

Which would run the simulation in a similar way and return simparams with an updated total_M to recreate the baseline prevalence.

Developer notes:

• Other paramterisation/calibration functions are currently in R/parameters.R.
• It would be appropriate to add a test in tests/testthat/test-biology.R to make sure that it indeed produces a pfpr2-10 that is close to the desired one.

[giovannic]

Note: if this relies on the ODE to calculate total_M, this may be impossible. Since adult mosquitoes are now modelled individually.

[pwinskill]

Hey @giovannic !

I've been playing around with a similar root-finding approach to the oneI implemented in the old model for calibrating a run. It is definitely not the most efficient approach, but is flexible, allowing a user to calibrate a specific model run to any number of outputpoints from any output variable. Would be good to hear if you think this kind of approach would be useful to implement in the package, or demonstrate in a vignette? It would be made more efficient if we can get a half decent starting guess for EIR, or the reduce the range of EIRs to search.

Here's a quick example:

library(malariasimulation)

# Objective (internal) function
objective <- function(eir, parameters, target_var, target, target_tt, summary_function, tolerance){
  message("Trying EIR: ", round(eir * 365, 4))
  p <- parameters %>%
    set_equilibrium(init_EIR = eir)
  # Setting seed - this might help??
  set.seed(1233445)
  raw <- malariasimulation::run_simulation(timesteps = timesteps, parameters = p)
  summary <- summary_function(raw)
  difference <- (summary[target_tt, target_var] - target)
  message("Current difference: ", round(difference, 3))
  difference[abs(difference) < tolerance] <- 0
  difference <- sum(difference)
  return(difference)
}

# Calibrate (user facing)
# target_var: Name of variable to calibrate to
# target: target values of target_var to calibrate to 
# target_tt: target timesteps at which to match target_var targets
# summary_function: a post-processing function, to allow second order outputs to be calibrated to
# tolerance: the tolerance require for calibration (average error on the same scale as target_var)
calibrate <- function(timesteps, parameters, target_var, target, target_tt, summary_function, tolerance){
  uniroot(objective,
          interval = c(0.01, 500) / 365,
          extendInt = "no",
          parameters = parameters,
          target_var = target_var,
          target = target,
          target_tt = target_tt,
          summary_function = summary_function,
          tolerance = tolerance, 
          trace = 1)
}


###  Example: ##################################################################
# Allow the user to define a processing function for second order output
sf <- function(x){
  x$prev_2_10 <- x$n_detect_730_3650 / x$n_730_3650
  return(x)
}
# Name the target variable to fit to
target_var <- "prev_2_10"
# Define target (in this case prev) to find
target <- c(0.1, 0.1)
# Time points at which to match target
target_tt <- c(500, 600)

p <- get_parameters(list(human_population = 5000))
timesteps <- 1000
out <- calibrate(timesteps = timesteps,
                 parameters = p,
                 target_var = target_var,
                 target = target,
                 target_tt = target_tt,
                 summary_function = sf,
                 tolerance = 0.001)
pars <- p %>%
  set_equilibrium(init_EIR = out$root)
raw <- malariasimulation::run_simulation(timesteps = timesteps, parameters = pars)
summary <- sf(raw)
plot(summary[,c("timestep", target_var)], t = "l", ylim = c(0, 0.3))
points(target ~ target_tt, col = "red", pch = 19)
################################################################################

[pwinskill]

Here's the output btw to that example. Fitting to two prevalence points (with reasonable tolerance):

[giovannic]

Closing, we would like to implement this in a separate package