Merge pull request #180 from mrc-ide/feat/time_vary_demog

Time varying demography

NAMESPACE

@@ -4,7 +4,6 @@ export(AL_params)
 export(DHA_PQP_params)
 export(SP_AQ_params)
 export(arab_params)
-export(find_birthrates)
 export(fun_params)
 export(gamb_params)
 export(get_correlation_parameters)

R/parameters.R

@@ -400,6 +400,8 @@ get_parameters <- function(overrides = list()) {
     severe_prevalence_rendering_max_ages = numeric(0),
     severe_incidence_rendering_min_ages = numeric(0),
     severe_incidence_rendering_max_ages = numeric(0),
+    age_group_rendering_min_ages = numeric(0),
+    age_group_rendering_max_ages = numeric(0),
     # misc
     custom_demography = FALSE,
     human_population = 100,

R/population_parameters.R

@@ -3,8 +3,6 @@
 #' @param parameters the model parameters
 #' @param agegroups vector of agegroups (in timesteps)
 #' @param timesteps vector of timesteps for each change in demography
-#' @param birthrates vector of new birthrates (number of individuals born per
-#' timestep) for each element of `timesteps`
 #' @param deathrates matrix of deathrates per age group per timestep.
 #' Rows are timesteps from the `timesteps` param. Columns are the age groups
 #' from the `agegroups` param.
@@ -13,49 +11,24 @@ set_demography <- function(
   parameters,
   agegroups,
   timesteps,
-  birthrates,
   deathrates
   ) {
-  stopifnot(all(agegroups > 0))
+
   stopifnot(all(timesteps >= 0))
-  stopifnot(all(birthrates > 0))
-  stopifnot(length(birthrates) == length(timesteps))
+  if(min(timesteps) != 0){
+    stop("Must include the baseline demography (timesteps must include 0),
+         when setting a custom demography")
+  }
+  stopifnot(all(agegroups > 0))
   stopifnot(all(deathrates > 0 & deathrates < 1))
   stopifnot(length(agegroups) == ncol(deathrates))
   stopifnot(length(timesteps) == nrow(deathrates))
   stopifnot(!is.unsorted(timesteps, strictly = TRUE))
-  stopifnot(length(timesteps) == 1) # changing population is not yet supported
+
   parameters$custom_demography <- TRUE
   parameters$deathrate_agegroups <- agegroups
   parameters$deathrate_timesteps <- timesteps
   parameters$deathrates <- deathrates
-  parameters$birthrate_timesteps <- timesteps
-  parameters$birthrates <- birthrates
-  n_age <- length(agegroups)
-
-  # set the populations
-  populations <- vapply(
-    seq(timesteps),
-    function(timestep) {
-      get_equilibrium_population(
-        agegroups,
-        birthrates[[timestep]],
-        parameters$deathrates[timestep,]
-      )
-    },
-    numeric(n_age)
-  )
-
-  deathrates <- vnapply(
-    seq(timesteps),
-    function(timestep) {
-      sum(
-         parameters$deathrates[timestep,] * populations[timestep,]
-      )
-    }
-  )
-  parameters$human_population <- round(colSums(populations))
-  parameters$human_population_timesteps <- timesteps
 
   parameters
 }
@@ -87,7 +60,6 @@ get_equilibrium_population <- function(age_high, birthrate, deathrates) {
 #' @param age_high a vector of age groups
 #' @param deathrates vector of deathrates for each age group
 #' @importFrom stats uniroot
-#' @export
 find_birthrates <- function(pops, age_high, deathrates) {
   vnapply(
     pops,
@@ -100,14 +72,6 @@ find_birthrates <- function(pops, age_high, deathrates) {
   )
 }
 
-get_birthrate <- function(parameters, timestep) {
-  if (!parameters$custom_demography) {
-    return(1 / parameters$average_age * get_human_population(parameters, timestep))
-  }
-  last_birthrate <- match_timestep(parameters$birthrate_timesteps, timestep)
-  parameters$birthrates[last_birthrate]
-}
-
 get_human_population <- function(parameters, timestep) {
   last_pop <- match_timestep(parameters$human_population_timesteps, timestep)
   parameters$human_population[last_pop]

R/processes.R

@@ -158,6 +158,11 @@ create_processes <- function(
       parameters,
       renderer
     ),
+    create_age_group_renderer(
+      variables$birth,
+      parameters,
+      renderer
+    ),
     create_compartmental_rendering_process(renderer, solvers, parameters)
   )
 

R/render.R

@@ -152,3 +152,22 @@ create_total_M_renderer_compartmental <- function(renderer, solvers, parameters)
     }
   }
 }
+
+create_age_group_renderer <- function(
+    birth,
+    parameters,
+    renderer
+) {
+  function(timestep) {
+    for (i in seq_along(parameters$age_group_rendering_min_ages)) {
+      lower <- parameters$age_group_rendering_min_ages[[i]]
+      upper <- parameters$age_group_rendering_max_ages[[i]]
+      in_age <- in_age_range(birth, timestep, lower, upper)
+      renderer$render(
+        paste0('n_age_', lower, '_', upper),
+        in_age$size(),
+        timestep
+      ) 
+    }
+  }
+}

R/utils.R

@@ -62,5 +62,6 @@ rtexp <- function(n, m, t) { itexp(runif(n), m, t) }
 #'@param t current timestep
 #'@noRd
 match_timestep <- function(ts, t) {
-  min(sum(ts < t) + 1, length(ts))
+  min(sum(ts <= t), length(ts))
 }
+

R/variables.R

@@ -406,11 +406,12 @@ calculate_initial_ages <- function(parameters) {
     )))
   }
 
+  deathrates <- parameters$deathrates[1, , drop = FALSE]
   age_high <- parameters$deathrate_agegroups
   age_width <- diff(c(0, age_high))
   age_low <- age_high - age_width
   n_age <- length(age_high)
-  birthrate <- get_birthrate(parameters, 0)
+  birthrate <- find_birthrates(parameters$human_population, age_high, deathrates)
   deathrates <- parameters$deathrates[1,]
 
   eq_pop <- get_equilibrium_population(age_high, birthrate, deathrates)

tests/testthat/test-demography.R

@@ -31,15 +31,14 @@ test_that('find_birthrates is consistent with get_equilibrium_population', {
 
 
 test_that('calculate_initial_ages calculates truncated exp custom demographic', {
-  parameters <- get_parameters()
+  parameters <- get_parameters(list(human_population = 4))
   ages <- c(50, 100) * 365
   deathrates <- c(.5, .75)
   parameters <- set_demography(
     parameters,
     agegroups = ages,
-    timesteps = 1,
-    birthrates = find_birthrates(4, ages, deathrates),
-    deathrates = matrix(deathrates, nrow=1, ncol=2)
+    timesteps = 0,
+    deathrates = matrix(deathrates, nrow = 1, ncol = 2)
   )
   mock_groups <- mockery::mock(c(2, 1, 2, 1))
   mock_rtexp <- mockery::mock(c(25 * 365, 30 * 365), c(25 * 365, 30 * 365))
@@ -56,3 +55,19 @@ test_that('calculate_initial_ages calculates truncated exp custom demographic',
   mockery::expect_args(mock_rtexp, 2, 2, .75, 50 * 365)
   expect_setequal(ages, c(25 * 365, 75 * 365, 30 * 365, 80 * 365))
 })
+
+test_that('match_timestep always gives 0 for constant demography', {
+  expect_equal(match_timestep(c(0), 0), 1)
+  expect_equal(match_timestep(c(0), 1), 1)
+  expect_equal(match_timestep(c(0), 50), 1)
+})
+
+test_that('match_timestep works on the boundaries', {
+  expect_equal(match_timestep(c(0, 50, 100), 0), 1)
+  expect_equal(match_timestep(c(0, 50, 100), 1), 1)
+  expect_equal(match_timestep(c(0, 50, 100), 49), 1)
+  expect_equal(match_timestep(c(0, 50, 100), 50), 2)
+  expect_equal(match_timestep(c(0, 50, 100), 99), 2)
+  expect_equal(match_timestep(c(0, 50, 100), 100), 3)
+  expect_equal(match_timestep(c(0, 50, 100), 101), 3)
+})

tests/testthat/test-mortality-integration.R

@@ -49,13 +49,13 @@ test_that('mortality_process resets humans correctly', {
 test_that('mortality_process samples deaths from a custom demography', {
   timestep <- 2
   parameters <- get_parameters()
+  parameters$human_population <- 4
   ages <- c(50, 100) * 365
   deaths <- c(.5, .75)
   parameters <- set_demography(
     parameters,
     agegroups = ages,
-    timesteps = 1,
-    birthrates = find_birthrates(4, ages, deaths),
+    timesteps = 0,
     deathrates = matrix(deaths, nrow=1, ncol=2)
   )
   events <- create_events(parameters)