From d0854db6fd3a135c043544eb479d7106b8e212c5 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Thu, 9 May 2024 10:39:22 +0100 Subject: [PATCH 1/2] Older age groups (older than equilibrium age object) are not correctly initiated. If outside of the age range, they get assigned 1 - the youngest age group. If we limit the max ages during these steps, we get something that is more accurate. --- R/variables.R | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/R/variables.R b/R/variables.R index 092a6431..6f6d9969 100644 --- a/R/variables.R +++ b/R/variables.R @@ -321,6 +321,8 @@ initial_immunity <- function( ) { if (!is.null(eq)) { age <- age / 365 + # older ages must be assigned within the equilibrium age bounds + age[age>=max(eq[[1]][, 'age'])] <- max(eq[[1]][, 'age'])-0.01 return(vnapply( seq_along(age), function(i) { @@ -338,6 +340,8 @@ initial_state <- function(parameters, age, groups, eq, states) { if (!is.null(eq)) { eq_states <- c('S', 'D', 'A', 'U', 'T') age <- age / 365 + # older ages must be assigned within the equilibrium age bounds + age[age>=max(eq[[1]][, 'age'])] <- max(eq[[1]][, 'age'])-0.01 return(vcapply( seq_along(age), function(i) { From 40e4f8785ff58364da21feb7e18c988a8e453d35 Mon Sep 17 00:00:00 2001 From: RJSheppard Date: Tue, 21 May 2024 12:52:45 +0100 Subject: [PATCH 2/2] Changed the initial age distribution to a truncated exponential to remove issues with state and immunity assignment when someone is older than the equilibrium solution. --- R/variables.R | 9 +++------ tests/testthat/test-demography.R | 9 +++++---- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/R/variables.R b/R/variables.R index 6f6d9969..86a13185 100644 --- a/R/variables.R +++ b/R/variables.R @@ -321,8 +321,6 @@ initial_immunity <- function( ) { if (!is.null(eq)) { age <- age / 365 - # older ages must be assigned within the equilibrium age bounds - age[age>=max(eq[[1]][, 'age'])] <- max(eq[[1]][, 'age'])-0.01 return(vnapply( seq_along(age), function(i) { @@ -340,8 +338,6 @@ initial_state <- function(parameters, age, groups, eq, states) { if (!is.null(eq)) { eq_states <- c('S', 'D', 'A', 'U', 'T') age <- age / 365 - # older ages must be assigned within the equilibrium age bounds - age[age>=max(eq[[1]][, 'age'])] <- max(eq[[1]][, 'age'])-0.01 return(vcapply( seq_along(age), function(i) { @@ -399,9 +395,10 @@ calculate_initial_ages <- function(parameters) { n_pop <- get_human_population(parameters, 0) # check if we've set up a custom demography if (!parameters$custom_demography) { - return(round(rexp( + return(round(rtexp( n_pop, - rate = 1 / parameters$average_age + 1 / parameters$average_age, + max(EQUILIBRIUM_AGES)*365 ))) } diff --git a/tests/testthat/test-demography.R b/tests/testthat/test-demography.R index 3b08dd37..ce760567 100644 --- a/tests/testthat/test-demography.R +++ b/tests/testthat/test-demography.R @@ -1,13 +1,14 @@ -test_that('calculate_initial_ages defaults to an exponential distribution', { +test_that('calculate_initial_ages defaults to a truncated exponential distribution', { parameters <- get_parameters(list(human_population = 4)) mock_exp <- mockery::mock(seq(4)) - mockery::stub(calculate_initial_ages, 'rexp', mock_exp) + mockery::stub(calculate_initial_ages, 'rtexp', mock_exp) ages <- calculate_initial_ages(parameters) mockery::expect_args( mock_exp, 1, parameters$human_population, - 1 / parameters$average_age + 1 / parameters$average_age, + max(EQUILIBRIUM_AGES)*365 ) }) @@ -70,4 +71,4 @@ test_that('match_timestep works on the boundaries', { 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) -}) \ No newline at end of file +})