Merge pull request #185 from mrc-ide/feat/pmc

PMC (IPTi) intervention

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: malariasimulation
 Title: An individual based model for malaria
-Version: 1.4.0
+Version: 1.4.1
 Authors@R: c(
   person(
     given = "Giovanni",

NAMESPACE

@@ -21,6 +21,7 @@ export(set_drugs)
 export(set_equilibrium)
 export(set_mass_rtss)
 export(set_mda)
+export(set_pmc)
 export(set_rtss_epi)
 export(set_smc)
 export(set_species)

R/correlation.R

@@ -1,4 +1,5 @@
 INTS <- c(
+  'pmc',
   'rtss',
   'mda',
   'smc',

R/mda_parameters.R

@@ -73,3 +73,33 @@ set_smc <- function(
   parameters$smc_max_ages <- max_ages
   parameters
 }
+
+#' @title Parameterise a perennial malaria chemoprevention (PMC, formerly IPIi)
+#' @param parameters a list of parameters to modify
+#' @param drug the index of the drug to administer
+#' @param timesteps a vector of timesteps for each round of PMC
+#' @param coverages a vector of the proportion of the target population who receive each
+#' round
+#' @param ages a vector of ages at which PMC is administered (in timesteps)
+#' @export
+set_pmc <- function(
+    parameters,
+    drug,
+    timesteps,
+    coverages,
+    ages
+) {
+
+  if(length(coverages) != length(timesteps)){
+    stop("coverages and timesteps do no align")
+  }
+  # check that the drug is valid
+  stopifnot((drug > 0) && (drug <= length(parameters$drug_rel_c)))
+
+  parameters$pmc <- TRUE
+  parameters$pmc_drug <- drug
+  parameters$pmc_timesteps <- timesteps
+  parameters$pmc_coverages <- coverages
+  parameters$pmc_ages <- ages
+  parameters
+}

R/parameters.R

@@ -182,7 +182,7 @@
 #' * rtss_ds - delay parameters for the antibody model, short-term weaning (mean and std. dev); default = 3.74502, 0.341185 (White MT et al. 2015 Lancet ID)
 #' * rtss_dl - delay parameters for the antibody model, long-term weaning (mean and std. dev); default = 6.30365, 0.396515 (White MT et al. 2015 Lancet ID)
 #'
-#' MDA and SMC parameters:
+#' MDA, SMC and PMC parameters:
 #' please set these parameters with the convenience functions in `mda_parameters.R`
 #'
 #' TBV parameters: 
@@ -372,6 +372,12 @@ get_parameters <- function(overrides = list()) {
     smc_coverages = NULL,
     smc_min_ages = -1,
     smc_max_ages = -1,
+    # PMC
+    pmc = FALSE,
+    pmc_drug = 0,
+    pmc_timesteps = NULL,
+    pmc_coverages = NULL,
+    pcs_ages = -1,
     # tbv
     tbv = FALSE,
     tbv_mt = 35,

R/pmc.R

@@ -0,0 +1,74 @@
+#' @title pmc process
+#'
+#' @description schedules individuals to be given perennial malaria 
+#' chemoprevention according to an age-based strategy
+#'
+#' @param variables list of variables in the model
+#' @param events a list of events in the model
+#' @param parameters the model parameters
+#' @param renderer a renderer object
+#' @param correlations correlation parameters
+#' @param coverages PMC coverage
+#' @param timesteps PMC coverage change timesteps
+#' @param drug PMC drug index
+#' @noRd
+create_pmc_process <- function(
+    variables,
+    events,
+    parameters,
+    renderer,
+    correlations,
+    coverages,
+    timesteps,
+    drug
+) {
+  function(timestep) {
+    timestep_index <- match_timestep(ts = timesteps, t = timestep)
+    if(timestep_index == 0){
+      return()
+    }
+    coverage <- coverages[timestep_index]
+    if(coverage == 0){
+      return()
+    }
+
+    age <- get_age(variables$birth$get_values(), timestep)
+
+    in_age <- which(age %in% parameters$pmc_ages)
+    target <- in_age[sample_intervention(in_age, 'pmc', coverage, correlations)]
+
+    successful_treatments <- bernoulli(
+      length(target),
+      parameters$drug_efficacy[[drug]]
+    )
+    to_move <- individual::Bitset$new(parameters$human_population)
+    to_move$insert(target[successful_treatments])
+    renderer$render('n_pmc_treated', to_move$size(), timestep)
+
+    if (to_move$size() > 0) {
+      # Move Diseased
+      diseased <- variables$state$get_index_of(c('D', 'A'))$and(to_move)
+      if (diseased$size() > 0) {
+        variables$state$queue_update('Tr', diseased)
+      }
+
+      # Move everyone else
+      other <- to_move$copy()$and(diseased$not(TRUE))
+      if (other$size() > 0) {
+        variables$state$queue_update('S', other)
+      }
+
+      # Update infectivity
+      variables$infectivity$queue_update(
+        variables$infectivity$get_values(
+          to_move
+        ) * parameters$drug_rel_c[[drug]],
+        to_move
+      )
+
+      # Update drug
+      variables$drug$queue_update(drug, to_move)
+      variables$drug_time$queue_update(timestep, to_move)
+    }
+  }
+}

R/processes.R

@@ -137,6 +137,25 @@ create_processes <- function(
       )
     )
   }
+
+  # =========
+  # PMC
+  # =========
+  if(!is.null(parameters$pmc_coverages)){
+    processes <- c(
+      processes,
+      create_pmc_process(
+        variables,
+        events,
+        parameters,
+        renderer,
+        correlations,
+        parameters$pmc_coverages,
+        parameters$pmc_timesteps,
+        parameters$pmc_drug
+      )
+    )
+  }
 
   # =========
   # Rendering

tests/testthat/test-pmc.R

@@ -0,0 +1,104 @@
+test_that("pmc parameterisation works", {
+  p <- get_parameters()
+
+  expect_false(p$pmc)
+  expect_equal(p$pmc_ages, NULL)
+  expect_equal(p$pmc_coverages, NULL)
+  expect_equal(p$pmc_timesteps, NULL)
+
+
+  p <- set_drugs(
+    parameters = p,
+    drugs = list(SP_AQ_params))
+  p <- set_pmc(
+    parameters = p,
+    drug = 1,
+    timesteps = c(5, 10),
+    coverages = c(0.5, 1),
+    ages = c(30, 60, 90)
+  )
+
+  expect_true(p$pmc)
+  expect_equal(p$pmc_ages, c(30, 60, 90))
+  expect_equal(p$pmc_coverages, c(0.5, 1))
+  expect_equal(p$pmc_timesteps, c(5, 10))
+  expect_equal(p$pmc_drug, 1)
+})
+
+
+test_that("pmc gives drugs to correct ages", {
+
+  p <- get_parameters(list(human_population = 6))
+  p <- set_drugs(
+    parameters = p,
+    drugs = list(SP_AQ_params))
+
+  p <- set_pmc(
+    parameters = p,
+    drug = 1,
+    timesteps = 10,
+    coverages = 1,
+    ages = c(30, 60, 90)
+  )
+  timestep <- 10
+  events <- create_events(p)
+  renderer <- individual:::Render$new(timestep)
+  variables <- create_variables(p)
+  variables$birth <- individual::IntegerVariable$new(
+    -c(10, 30, 60, 89, 90, 36500) + 10
+  )
+  variables$state <- mock_category(
+    c('D', 'S', 'A', 'U', 'Tr'),
+    c('D', 'S', 'A', 'U', 'D', 'S')
+  )
+  variables$drug <- mock_integer(rep(0, 6))
+  variables$drug_time <- mock_integer(rep(-1, 6))
+  mockery::stub(sample_intervention, 'bernoulli', mockery::mock(c(TRUE, TRUE, TRUE)))
+
+  process <- create_pmc_process(
+    variables = variables,
+    events = events, 
+    parameters = p,
+    renderer = renderer,
+    correlations = get_correlation_parameters(p),
+    coverages = p$pmc_coverages,
+    timesteps = p$pmc_timesteps,
+    drug = p$pmc_drug
+  )
+  # mock the treatment success
+  mockery::stub(process, 'bernoulli', mockery::mock(c(TRUE, TRUE, TRUE)))
+  process(timestep)
+
+  # Three treatments given
+  expect_equal(renderer$to_dataframe(),
+               data.frame(timestep = 1:10,
+                          n_pmc_treated = c(rep(NA, 9), 3)))
+
+  # Individuals 3 and 5, are correct age and in D or A states
+  expect_bitset_update(
+    variables$state$queue_update_mock(),
+    'Tr',
+    c(3, 5),
+    1
+  )
+  # Individual 2 is correct age and in S state
+  expect_bitset_update(
+    variables$state$queue_update_mock(),
+    'S',
+    2,
+    2
+  )
+  # Drug is recorded as given
+  expect_bitset_update(
+    variables$drug$queue_update_mock(),
+    1,
+    c(2, 3, 5)
+  )
+  # Drug time is recorded
+  expect_bitset_update(
+    variables$drug_time$queue_update_mock(),
+    10,
+    c(2, 3, 5)
+  )
+})
+