inconsistency in set_drug vector lengths does not introduce error

[htopazian]

I have identified an area where a small test could be introduced to alert the user if the parameters involved in set_drugs are not consistent in length (in this case drug_prophylaxis_shape and drug_prophylaxis_scale do not match drug_efficacyand drug_rel_c). This problem can lead to incorrect results where the model is implementing a different set of drug characteristics than the user thinks they are setting.

I have included an example below! Perhaps we could put in a validation test to error out if these vector lengths are inconsistent?

library(tidyverse)
library(malariasimulation)


# set parameters
year <- 365
month <- year / 12
warmup <- 3 * year
sim_length <- 5 * year


# starting parameters
params <- get_parameters(list(
  human_population = 2000,
  model_seasonality = FALSE,
  individual_mosquitoes = FALSE))

# progress bar
params$progress_bar <- TRUE

# outputs
params$clinical_incidence_rendering_min_ages = c(0)
params$clinical_incidence_rendering_max_ages = c(100 * year)

# drugs
params <- set_drugs(
  parameters = params,
  list(AL_params, SP_AQ_params, DHA_PQP_params))

# glance at shape and scale
params$drug_prophylaxis_scale; params$drug_prophylaxis_shape

# replace shape and scale with a vector of 2 (not 3 as it should be)
# AL default, SP: https://doi.org/10.1016/S2214-109X(22)00416-8 supplement
params$drug_prophylaxis_scale <- c(10.6, 39.34) 
params$drug_prophylaxis_shape <- c(11.3, 3.40) 

# glance at shape and scale
params$drug_prophylaxis_scale; params$drug_prophylaxis_shape

# MDA
first <- warmup + 100 # first dose
params <- set_mda(
  parameters = params,
  drug = 3, # DHA_PQP for MDA 
  timesteps = c(first, first + round(1 * month), first + round(2 * month)),
  coverages = rep(0.80, 3), 
  min_ages = rep(0, 3), 
  max_ages = rep(100 * year, 3)
)

# set equilibrium
params <- set_equilibrium(params, 100)

# run simulation
set.seed(123)

output <- run_simulation(
  timesteps = warmup + sim_length,
  parameters = params) |>
  mutate(timestep = timestep - warmup,
         year = ceiling(timestep / year)) |>
  group_by(year) |>
  summarize(n_inc_clinical_0_36500 = sum(n_inc_clinical_0_36500),
            n_age_0_36500 = mean(n_age_0_36500)) |>
  mutate(clinical_incidence = n_inc_clinical_0_36500 / n_age_0_36500,
         run = "wrong vector length")


# run second model with drug vector fix
params$drug_prophylaxis_scale <- c(10.6, 39.34, 28.1) 
params$drug_prophylaxis_shape <- c(11.3, 3.40, 4.4) 

set.seed(123)

output2 <- run_simulation(
  timesteps = warmup + sim_length,
  parameters = params) |>
  mutate(timestep = timestep - warmup,
         year = ceiling(timestep / year)) |>
  group_by(year) |>
  summarize(n_inc_clinical_0_36500 = sum(n_inc_clinical_0_36500),
            n_age_0_36500 = mean(n_age_0_36500)) |>
  mutate(clinical_incidence = n_inc_clinical_0_36500 / n_age_0_36500,
         run = "correct vector length")

# plot
ggplot(data = rbind(output, output2)) +
  geom_line(aes(x = year, y = clinical_incidence, color = run))