Function for setting age-structured outputs for incidence, prevalence…

… and immunity. Discrete age groups are set using a list, where vectors within the list result in contiguous age group rendering. This is a correction attempting to remove issues around line endings.
mrc-ide · May 8, 2024 · 80eb376 · 80eb376
1 parent 1111500
commit 80eb376
Show file tree

Hide file tree

Showing 5 changed files with 156 additions and 35 deletions.
diff --git a/R/output_parameters.R b/R/output_parameters.R
@@ -0,0 +1,39 @@
+#' @title Parameterise age grouped output rendering
+#'
+#' @details this function produces discrete and contiguous age groups, inclusive of the lower
+#' age limit and exclusive of the upper age limit: e.g., list(c(0, 10, 100), c(200, 250) will produce
+#' three age groups: 0-9, 10-99 and 200-249 in days.
+#' @param parameters the model parameters
+#' @param age_group age breaks for population size outputs; default = NULL
+#' @param incidence age breaks for incidence outputs (D+Tr+A); default = NULL
+#' @param clinical_incidence age breaks for clinical incidence outputs (symptomatic); default = c(0, 1825)
+#' @param severe_incidence age breaks for severe incidence outputs; default = NULL
+#' @param prevalence age breaks for clinical prevalence outputs (pcr and lm detectable infections); default = c(730, 3650)
+#' @export
+#'
+set_epi_outputs <- function(parameters,
+                            age_group = NULL,
+                            incidence = NULL,
+                            clinical_incidence = NULL,
+                            severe_incidence = NULL,
+                            prevalence = NULL,
+                            ica = NULL,
+                            icm = NULL,
+                            iva = NULL,
+                            ivm = NULL,
+                            id = NULL,
+                            ib = NULL
+){
+
+  parent_formals <- names(formals())
+  parent_formals <- parent_formals[which(parent_formals != "parameters")]
+  outputs <- parent_formals[!unlist(lapply(parent_formals, function(x){is.null(get(x))}))]
+
+  for (output in outputs) {
+    if(!is.list(get(output))){stop("Each age input must be a list of vectors")}
+    parameters[[paste0(output, "_rendering_min_ages")]] <- unlist(lapply(get(output), function(x){x[-length(x)]}))
+    parameters[[paste0(output, "_rendering_max_ages")]] <- unlist(lapply(get(output), function(x){x[-1]-1}))
+  }
+
+  parameters
+}
diff --git a/R/processes.R b/R/processes.R
@@ -30,7 +30,7 @@ create_processes <- function(
     correlations,
     lagged_eir,
     lagged_infectivity,
-    timesteps, 
+    timesteps,
     mixing = 1,
     mixing_index = 1
 ) {
@@ -48,7 +48,7 @@ create_processes <- function(
     create_exponential_decay_process(variables$iva, parameters$rva),
     create_exponential_decay_process(variables$id, parameters$rid)
   )
-  
+
   if (parameters$individual_mosquitoes) {
     processes <- c(
       processes,
@@ -61,7 +61,7 @@ create_processes <- function(
       )
     )
   }
-  
+
   # ==============================
   # Biting and mortality processes
   # ==============================
@@ -106,22 +106,22 @@ create_processes <- function(
       0
     )
   )
-  
+
   # =======================
   # Antimalarial Resistance
   # =======================
   # Add an a new process which governs the transition from Tr to S when
   # antimalarial resistance is simulated. The rate of transition switches
   # from a parameter to a variable when antimalarial resistance == TRUE.
-  
+
   # Assign the dt input to a separate object with the default single parameter value:
   dt_input <- parameters$dt
-  
-  # If antimalarial resistance is switched on, assign dt variable values to the 
+
+  # If antimalarial resistance is switched on, assign dt variable values to the
   if(parameters$antimalarial_resistance) {
     dt_input <- variables$dt
   }
-  
+
   # Create the progression process for Tr --> S specifying dt_input as the rate:
   processes <- c(
     processes,
@@ -134,15 +134,15 @@ create_processes <- function(
       0
     )
   )
-  
+
   # ===============
   # ODE integration
   # ===============
   processes <- c(
     processes,
     create_solver_stepping_process(solvers, parameters)
   )
-  
+
   # =========
   # RTS,S EPI
   # =========
@@ -159,7 +159,7 @@ create_processes <- function(
       )
     )
   }
-  
+
   # =========
   # PMC
   # =========
@@ -178,10 +178,13 @@ create_processes <- function(
       )
     )
   }
-  
+
   # =========
   # Rendering
   # =========
+
+  imm_var_names <- c('ica', 'icm', 'id', 'ib', 'iva', 'ivm')
+
   processes <- c(
     processes,
     individual::categorical_count_renderer_process(
@@ -191,10 +194,10 @@ create_processes <- function(
     ),
     create_variable_mean_renderer_process(
       renderer,
-      c('ica', 'icm', 'ib', 'id', 'iva', 'ivm'),
-      variables[c('ica', 'icm', 'ib', 'id', 'iva', 'ivm')]
+      imm_var_names,
+      variables[imm_var_names]
     ),
-    create_prevelance_renderer(
+    create_prevalence_renderer(
       variables$state,
       variables$birth,
       variables$id,
@@ -206,9 +209,16 @@ create_processes <- function(
       parameters,
       renderer
     ),
+    create_age_variable_mean_renderer_process(
+      imm_var_names[paste0(imm_var_names,"_rendering_min_ages") %in% names(parameters)],
+      variables[imm_var_names[paste0(imm_var_names,"_rendering_min_ages") %in% names(parameters)]],
+      variables$birth,
+      parameters,
+      renderer
+    ),
     create_compartmental_rendering_process(renderer, solvers, parameters)
   )
-  
+
   if (parameters$individual_mosquitoes) {
     processes <- c(
       processes,
@@ -230,11 +240,11 @@ create_processes <- function(
       )
     )
   }
-  
+
   # ======================
   # Intervention processes
   # ======================
-  
+
   if (parameters$bednets) {
     processes <- c(
       processes,
@@ -247,14 +257,14 @@ create_processes <- function(
       net_usage_renderer(variables$net_time, renderer)
     )
   }
-  
+
   if (parameters$spraying) {
     processes <- c(
       processes,
       indoor_spraying(variables$spray_time, parameters, correlations)
     )
   }
-  
+
   # ======================
   # Progress bar process
   # ======================
@@ -264,7 +274,7 @@ create_processes <- function(
       create_progress_process(timesteps)
     )
   }
-  
+
   processes
 }
 
@@ -289,7 +299,7 @@ create_exponential_decay_process <- function(variable, rate) {
 #' @title Create and initialise lagged_infectivity object
 #'
 #' @param variables model variables for initialisation
-#' @param parameters model parameters 
+#' @param parameters model parameters
 #' @noRd
 create_lagged_infectivity <- function(variables, parameters) {
   age <- get_age(variables$birth$get_values(), 0)
@@ -306,7 +316,7 @@ create_lagged_infectivity <- function(variables, parameters) {
 #'
 #' @param variables model variables for initialisation
 #' @param solvers model differential equation solvers
-#' @param parameters model parameters 
+#' @param parameters model parameters
 #' @noRd
 create_lagged_eir <- function(variables, solvers, parameters) {
   lapply(

diff --git a/R/render.R b/R/render.R
@@ -3,18 +3,18 @@ in_age_range <- function(birth, timestep, lower, upper) {
 }
 
 #' @title Render prevalence statistics
-#' 
+#'
 #' @description renders prevalence numerators and denominators for indivduals
 #' detected by microscopy and with severe malaria
-#' 
+#'
 #' @param state human infection state
 #' @param birth variable for birth of the individual
 #' @param immunity to detection
 #' @param parameters model parameters
 #' @param renderer model renderer
-#' 
+#'
 #' @noRd
-create_prevelance_renderer <- function(
+create_prevalence_renderer <- function(
   state,
   birth,
   immunity,
@@ -41,7 +41,7 @@ create_prevelance_renderer <- function(
         paste0('n_', lower, '_', upper),
         in_age$size(),
         timestep
-      ) 
+      )
       renderer$render(
         paste0('n_detect_', lower, '_', upper),
         in_age$copy()$and(detected)$size(),
@@ -59,9 +59,9 @@ create_prevelance_renderer <- function(
 }
 
 #' @title Render incidence statistics
-#' 
+#'
 #' @description renders incidence (new for this timestep) for indivduals
-#' 
+#'
 #' @param birth variable for birth of the individual
 #' @param renderer object for model outputs
 #' @param target incidence population
@@ -71,7 +71,7 @@ create_prevelance_renderer <- function(
 #' @param lowers age bounds
 #' @param uppers age bounds
 #' @param timestep current target
-#' 
+#'
 #' @noRd
 incidence_renderer <- function(
   birth,
@@ -119,6 +119,29 @@ create_variable_mean_renderer_process <- function(
   }
 }
 
+create_age_variable_mean_renderer_process <- function(
+    names,
+    variables,
+    birth,
+    parameters,
+    renderer
+) {
+  function(timestep) {
+    for (i in seq_along(variables)) {
+      for (j in seq_along(parameters[[paste0(names[[i]],"_rendering_min_ages")]])) {
+        lower <- parameters[[paste0(names[[i]],"_rendering_min_ages")]][[j]]
+        upper <- parameters[[paste0(names[[i]],"_rendering_max_ages")]][[j]]
+        in_age <- in_age_range(birth, timestep, lower, upper)
+        renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+        renderer$render(
+          paste0(names[[i]], '_mean_', lower, '_', upper),
+          mean(variables[[i]]$get_values(index = in_age)),
+          timestep
+        )
+      }
+    }
+  }
+}
 create_vector_count_renderer_individual <- function(
   mosquito_state,
   species,
@@ -168,7 +191,7 @@ create_age_group_renderer <- function(
         paste0('n_age_', lower, '_', upper),
         in_age$size(),
         timestep
-      ) 
+      )
     }
   }
 }
diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R
@@ -11,7 +11,7 @@ test_that('that default rendering works', {
   immunity <- individual::DoubleVariable$new(rep(1, 4))
 
   renderer <- mock_render(1)
-  process <- create_prevelance_renderer(
+  process <- create_prevalence_renderer(
     state,
     birth,
     immunity,
@@ -30,7 +30,7 @@ test_that('that default rendering works', {
     3,
     timestep
   )
-  
+
   mockery::expect_args(
     renderer$render_mock(),
     2,
@@ -62,7 +62,7 @@ test_that('that default rendering works when no one is in the age range', {
   immunity <- individual::DoubleVariable$new(rep(1, 4))
 
   renderer <- mock_render(1)
-  process <- create_prevelance_renderer(
+  process <- create_prevalence_renderer(
     state,
     birth,
     immunity,

diff --git a/tests/testthat/test_output.R b/tests/testthat/test_output.R
@@ -0,0 +1,49 @@
+test_that('Test age parameter function works', {
+  parameters <- get_parameters()
+  age_limits <- c(0,1,2,3)*365
+  parameters <- set_epi_outputs(parameters,
+                                age_group = list(age_limits),
+                                incidence = list(age_limits+1),
+                                clinical_incidence = list(age_limits+2),
+                                severe_incidence = list(age_limits+3),
+                                prevalence = list(age_limits+4),
+                                ica = list(age_limits+5),
+                                icm = list(age_limits+6),
+                                id = list(age_limits+7),
+                                ib = list(age_limits+8),
+                                iva = list(age_limits+9),
+                                ivm = list(age_limits+10)
+  )
+
+  sim <- run_simulation(timesteps = 1, parameters)
+
+  expect_true(
+    all(
+      paste0(rep(c("n_age",
+                   "n", "n_inc", "p_inc",
+                   "n","n_inc_clinical","p_inc_clinical",
+                   "n","n_inc_severe","p_inc_severe",
+                   "n","n_detect","p_detect",
+                   "ica_mean", "icm_mean","id_mean","ib_mean","iva_mean","ivm_mean"), each = 3),"_",
+             age_limits[-4]+rep(c(0,rep(c(1:4), each = 3),5:10), each = 3),"_",age_limits[-1]-1+rep(c(0,rep(c(1:4), each = 3),5:10), each = 3)
+      ) %in%
+        names(sim)))
+
+})
+
+test_that('Test age parameter multiple sequences function works', {
+  parameters <- get_parameters()
+  age_limits <- c(0,1,2,3)*365
+  parameters <- set_epi_outputs(parameters,
+                                age_group = list(age_limits, age_limits+100))
+
+  sim <- run_simulation(timesteps = 1, parameters)
+
+  expect_true(
+    all(
+      paste0(rep(c("n_age"), each = 3),"_",
+             c(age_limits[-4], age_limits[-4]+100), "_", c(age_limits[-1], age_limits[-1]+100)-1
+      ) %in%
+        names(sim)))
+
+})