From 25fa2aa81b828baf9893040a600eacc1f52d9a4d Mon Sep 17 00:00:00 2001 From: ruthmccabe Date: Sun, 6 Oct 2024 15:56:39 +0100 Subject: [PATCH 1/7] vaccine effectiveness in as default param given now explicitly modelling 1 and 2 doses --- R/parameters.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/R/parameters.R b/R/parameters.R index ed46d77..db41f5f 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -238,7 +238,13 @@ parameters_fixed <- function(region, initial_infections, overrides = list()) { }else{ stop("combination of population size and seeding infections is incompatible, please review.") } - + + ## vaccine efficacy + ve_I <- c(0.736,0,0.736,0.818) + if(length(ve_I)!=n_vax){ + stop("incorrect dimension for ve_I") + } + params_list = list( region = region, n_group = n_group, @@ -266,7 +272,7 @@ parameters_fixed <- function(region, initial_infections, overrides = list()) { # vaccination_coverage_target = matrix(0.01, nrow = n_group, ncol = N_prioritisation_steps), vaccine_uptake = rep(0.8, n_group), ve_T = rep(0, n_vax), - ve_I = rep(0, n_vax), + ve_I = ve_I, vaccination_coverage_target_1st_dose_prop = 0.8, vaccination_coverage_target_2nd_dose_prop = 0.5, vaccination_campaign_length = vaccination_campaign_length, From b12dbc2aec1cd6c68cfc1a80daa159c458aed835 Mon Sep 17 00:00:00 2001 From: ruthmccabe Date: Sun, 6 Oct 2024 16:04:23 +0100 Subject: [PATCH 2/7] vaccine-stratified CFR @lwhittles need you to have a look --- R/parameters.R | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/R/parameters.R b/R/parameters.R index db41f5f..530fed2 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -217,12 +217,20 @@ parameters_fixed <- function(region, initial_infections, overrides = list()) { # CFR from Whittles 2024, 5-year bands to 40 - CFR <- rep(0, n_group) - names(CFR) <- names(demographic_params$N0) - CFR[which(age_bins$end < 40)] <- c(0.102, 0.054, 0.035, 0.026, 0.02, 0.016, 0.013, 0.012) - CFR[which(age_bins$start >= 40)] <- 0.01 - CFR["SW"] <- CFR["20-24"] - CFR["PBS"] <- CFR["35-39"] + CFR <- matrix(0, nrow = n_group, ncol = n_vax) + rownames(CFR) <- names(demographic_params$N0) + + ## unvaccinated + CFR[which(age_bins$end < 40),2] <- c(0.102, 0.054, 0.035, 0.026, 0.02, 0.016, 0.013, 0.012) + CFR[which(age_bins$start >= 40),2] <- 0.01 + CFR["SW",2] <- CFR["20-24",2] + CFR["PBS",2] <- CFR["35-39",2] + + ## vaccinated (currently no distinction, LW to review) + CFR[which(age_bins$end < 40),c(1,3,4)] <- c(0.048, 0.025, 0.016, 0.012, 0.009, 0.007, 0.006, 0.005) + CFR[which(age_bins$start >= 40),c(1,3,4)] <- 0.004 + CFR["SW",c(1,3,4)] <- CFR["20-24",c(1,3,4)] + CFR["PBS",c(1,3,4)] <- CFR["35-39",c(1,3,4)] vaccination_campaign_length <- 1 @@ -265,7 +273,7 @@ parameters_fixed <- function(region, initial_infections, overrides = list()) { gamma_E = 1 / 7, # 1/7 based on Besombes et al. on 29 clade I patients gamma_Ir = 1 / 18, # Jezek 1988 "clinical features of 282.." gamma_Id = 1 / 10, # Jezek 1988 - CFR = matrix(CFR, nrow = n_group, ncol = n_vax, byrow = FALSE), + CFR = as.matrix(CFR), m_sex = demographic_params$m_sex, m_gen_pop = demographic_params$m_gen_pop, prioritisation_strategy = matrix(1, nrow = n_group, ncol = N_prioritisation_steps), From 0a2baa717239d9e029fe6a0f04aece3fa4733a88 Mon Sep 17 00:00:00 2001 From: ruthmccabe Date: Sun, 6 Oct 2024 16:08:44 +0100 Subject: [PATCH 3/7] bumped version number --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 13d38c8..b9b27cf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: mpoxseir Title: Stochastic compartmental model of mpox transmission -Version: 0.1.0 +Version: 0.1.1 Authors@R: c(person("Lilith", "Whittles", role = c("aut", "cre"), email = "l.whittles@imperial.ac.uk"), person("Ruth", "McCabe", role = c("aut")), From 0ffe4e4a0c73e6ec92e36d699c767ca354d4aa11 Mon Sep 17 00:00:00 2001 From: Lilith Whittles Date: Fri, 11 Oct 2024 16:30:59 +0100 Subject: [PATCH 4/7] set vaccine parameters to 0 in default params - we will vary this in onward simulation add helper functions to get compartment indices pull out use_ve_D as a parameter in parameters_fixed so that we can do fits allowing for different CFR by vaccine status, or not as required --- R/parameters.R | 166 +++++++++++++++++++++---------- R/run.R | 29 +++--- tests/testthat/helper-mpoxseir.R | 88 ++++------------ 3 files changed, 148 insertions(+), 135 deletions(-) diff --git a/R/parameters.R b/R/parameters.R index 530fed2..55e3bbf 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -70,10 +70,11 @@ parameters_demographic <- function() { marginalise <- function(M) (rowSums(M) + diag(M)) / 2 + compartment_indices <- get_compartment_indices() n_age <- length(N_age) idx_age <- seq_len(n_age) - n_group <- n_age + 2 - nms_group <- c(age_bins$label, "SW", "PBS") + n_group <- compartment_indices$dim$group + nms_group <- names(compartment_indices$group) ### General population @@ -99,8 +100,12 @@ parameters_demographic <- function() { ## set up sexual contact matrix for parameterisation in transform function m_sex <- matrix(0, n_group, n_group, dimnames = list(nms_group, nms_group)) - # proportion of susceptibles estimated to have smallpox vaccine - sus_prop <- c(rep(1,8),0.54,0.29,0.29,0.23,0.21,0.21,0.21,0.21,1,1) + # proportion of susceptibles estimated to be unvaccinated (historically) + p_unvaccinated <- setNames(rep(0, n_group), nms_group) + p_unvaccinated[which(age_bins$end < 40)] <- 1 + p_unvaccinated[which(age_bins$start >= 40)] <- + c(0.54, 0.29, 0.29, 0.23, 0.21, 0.21, 0.21, 0.21) + p_unvaccinated["SW"] <- p_unvaccinated["PBS"] <- 1 # province populations province_pop = list("equateur" = 1712000, @@ -114,8 +119,8 @@ parameters_demographic <- function() { m_sex = m_sex, total_contacts_gen_pop = M_all, #total_contacts_sex = M_sex, - n_vax = 4, - sus_prop = sus_prop, + n_vax = compartment_indices$dim$vax, + p_unvaccinated = p_unvaccinated, province_pop = province_pop ) } @@ -140,6 +145,39 @@ get_age_bins <- function() { data.frame(label = label, start = start, end = end) } + + +##' A function that gets the compartment indices used in the model +##' +##' @title Get compartment indices used throughout the model +##' +##' @return A list containing entries: `dim` giving the dimensions of each +##' compartment, currently: group = 18, vax = 4; +##' `group`, a named list of the array indices corresponding to the +##' first dimension of model compartments: 16 5-year age bands + SW + PBS; +##' and `vax`, a named list of the vaccine strata used in +##' the second dimension of model compartments: +##' 1. historic smallpox; 2. unvaccinated; 3. one-dose; 4. two-dose. +##' +##' Use this function wherever you need to refer to these standards in your code +##' to avoid duplication, and make modifying universally easier. +##' +##' @export +get_compartment_indices <- function() { + age_bins <- get_age_bins() + n_group <- nrow(age_bins) + 2 + groups <- seq_len(n_group) + names(groups) <- c(age_bins$label, "SW", "PBS") + + n_vax <- 4 + vax_strata <- seq_len(n_vax) + names(vax_strata) <- c("historic", "unvaccinated", "one_dose", "two_dose") + + list(dim = list(group = n_group, vax = n_vax), + group = as.list(groups), + vax = as.list(vax_strata)) +} + # Function to allocate N individuals into m groups based on weights assign_seeds <- function(N, w) { @@ -167,7 +205,8 @@ assign_seeds <- function(N, w) { ##' `"sudkivu"` ##' ##' @param initial_infections The initial number of infections -##' +##' @param use_ve_D logical, indicating whether model should allow for vaccine +##' efficacy against death (above and beyond protection against infection) ##' @param overrides A list, containing any parameters for which you want to ##' override the default values ##' @@ -176,7 +215,7 @@ assign_seeds <- function(N, w) { ##' @export ##' #' @export -parameters_fixed <- function(region, initial_infections, overrides = list()) { +parameters_fixed <- function(region, initial_infections, use_ve_D = FALSE, overrides = list()) { ## Checking region if (!(region %in% c("equateur", "sudkivu"))) { @@ -186,72 +225,95 @@ parameters_fixed <- function(region, initial_infections, overrides = list()) { ## Initialising variable that other parameters depend on demographic_params <- parameters_demographic() age_bins <- get_age_bins() + idx_compartment <- get_compartment_indices() n_group <- demographic_params$n_group n_vax <- demographic_params$n_vax + idx_unvax <- idx_compartment$vax$unvaccinated + idx_historic_vax <- idx_compartment$vax$historic + + ## standard vaccination parameters, we may want to move this entire section to + ## the transform in future if we want to allow for vaccine uncertainty + ## VE against infection + ve_I <- c(0.736, 0, 0.736, 0.818) # Berry et al. + ## VE against onward transmission + ve_T <- rep(0, n_vax) + + N <- demographic_params$province_pop[[region]] N0 <- round(N * demographic_params$N0 / sum(demographic_params$N0)) # total number in each age-group N_prioritisation_steps <- 1 + RR_z <- c(0.977, 1, 0.444, rep(0.078, n_group - 3)) # Jezek 1988 zoonotic + Jezek 1987 + ## Seed infections in the unvaccinated group in a region-specific manner X0 <- matrix(0, nrow = n_group, ncol = n_vax) - Ea0 <- matrix(0, nrow = n_group, ncol = n_vax) - RR_z <- c(0.977, 1, 0.444, rep(0.078, n_group - 3)) # Jezek 1988 zoonotic + Jezek 1987 + Ea0 <- X0 + + if (region == "sudkivu") { # seeding in sex workers in Sud Kivu ## Extract sex-worker index and put initial infections in this group (unvaccinated strata) - sw_index <- which(colnames(demographic_params$m_gen_pop) == "SW") - Ea0[sw_index, 2] <- initial_infections + index_sw <- which(colnames(demographic_params$m_gen_pop) == "SW") + Ea0[index_sw, idx_unvax] <- initial_infections } else if (region == "equateur") { # seeding in general pop in proportion to zoonotic risk in equateur ## Extract gen-pop index and put initial infections in this group (unvaccinated strata) in proportion to zoonotic risk - gen_pop_index <- seq_len(nrow(age_bins)) - seeding_infections <- assign_seeds(initial_infections, w = RR_z[gen_pop_index]) + index_gen_pop <- seq_len(nrow(age_bins)) + seeding_infections <- assign_seeds(initial_infections, w = RR_z[index_gen_pop]) - Ea0[gen_pop_index, 2] <- seeding_infections + Ea0[index_gen_pop, idx_unvax] <- seeding_infections } else { stop("something is wrong with the name of the region - change to sudkivu or equateur") } + + p_unvaccinated <- demographic_params$p_unvaccinated + ## update with assignment into smallpox vaccination compartment (j = 1) + ## if we have a small population then seeding with 5 infections per compt and + ## strata could be a problem (e.g. with N = 10,000 too few SW) + ## split population into historically- / un-vaccinated first before seeding + ## need to think about how this will interact with equilibrium position for + ## seeding in endemic areas. + S0 <- X0 + S0[, idx_unvax] <- round(N0 * p_unvaccinated) + S0[, idx_historic_vax] <- N0 - S0[, idx_unvax] + S0 <- S0 - Ea0 + if (any (S0 < 0)) { + stop("population size and seeding infections is incompatible") + } + CFR_unvax <- rep(0, n_group) + names(CFR_unvax) <- names(demographic_params$N0) # CFR from Whittles 2024, 5-year bands to 40 - - CFR <- matrix(0, nrow = n_group, ncol = n_vax) - rownames(CFR) <- names(demographic_params$N0) - - ## unvaccinated - CFR[which(age_bins$end < 40),2] <- c(0.102, 0.054, 0.035, 0.026, 0.02, 0.016, 0.013, 0.012) - CFR[which(age_bins$start >= 40),2] <- 0.01 - CFR["SW",2] <- CFR["20-24",2] - CFR["PBS",2] <- CFR["35-39",2] + CFR_unvax[which(age_bins$end < 40)] <- + c(0.102, 0.054, 0.035, 0.026, 0.02, 0.016, 0.013, 0.012) + CFR_unvax[which(age_bins$start >= 40)] <- 0.01 + + # Create matrix [group | vax strata] + # Vax strata: 1. historic smallpox; 2. unvaccinated; 3. one-dose; 4. two-dose + # initially all strata contain CFR by age for unvaccinated, which we then + # modify to incorporate any VE against death + CFR <- matrix(CFR_unvax, nrow = n_group, ncol = n_vax) + rownames(CFR) <- names(CFR_unvax) - ## vaccinated (currently no distinction, LW to review) - CFR[which(age_bins$end < 40),c(1,3,4)] <- c(0.048, 0.025, 0.016, 0.012, 0.009, 0.007, 0.006, 0.005) - CFR[which(age_bins$start >= 40),c(1,3,4)] <- 0.004 - CFR["SW",c(1,3,4)] <- CFR["20-24",c(1,3,4)] - CFR["PBS",c(1,3,4)] <- CFR["35-39",c(1,3,4)] + if (use_ve_D) { + # If allowing for vaccine protection against death + CFR_historic_vax <- CFR_unvax + CFR_historic_vax[which(age_bins$end < 40)] <- + c(0.048, 0.025, 0.016, 0.012, 0.009, 0.007, 0.006, 0.005) + CFR_historic_vax[which(age_bins$start >= 40)] <- 0.004 + + # use same CFR for all vaccinated regardless of efficacy + CFR[, -idx_unvax] <- CFR_historic_vax + } - vaccination_campaign_length <- 1 + CFR["SW", ] <- CFR["20-24", ] + CFR["PBS", ] <- CFR["35-39", ] - ## update with assignment into smallpox vaccination compartment (j=1) - ## if we have a small population then seeding with 5 infections per compt and strata could be a problem (e.g. with N=10000 too few SW) - if(all(rowSums(Ea0) Date: Fri, 11 Oct 2024 16:31:07 +0100 Subject: [PATCH 5/7] update docs --- NAMESPACE | 1 + man/create_transform_params.Rd | 14 +++++++++++--- man/get_compartment_indices.Rd | 25 +++++++++++++++++++++++++ man/parameters_fixed.Rd | 10 +++++++++- 4 files changed, 46 insertions(+), 4 deletions(-) create mode 100644 man/get_compartment_indices.Rd diff --git a/NAMESPACE b/NAMESPACE index d58f68e..ed36dd3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,6 +2,7 @@ export(create_transform_params) export(get_age_bins) +export(get_compartment_indices) export(model_index) export(model_targeted_vax) export(parameters_demographic) diff --git a/man/create_transform_params.Rd b/man/create_transform_params.Rd index d22359d..f7376b8 100644 --- a/man/create_transform_params.Rd +++ b/man/create_transform_params.Rd @@ -4,16 +4,24 @@ \alias{create_transform_params} \title{Create transform function} \usage{ -create_transform_params(region, initial_infections, overrides = list()) +create_transform_params( + region, + initial_infections, + use_ve_D = FALSE, + overrides = list() +) } \arguments{ -\item{region}{The region for fitting, must be either \code{"equateur"} or +\item{region}{The region for the parameters, must be either \code{"equateur"} or \code{"sudkivu"}} \item{initial_infections}{The initial number of infections} +\item{use_ve_D}{logical, indicating whether model should allow for vaccine +efficacy against death (above and beyond protection against infection)} + \item{overrides}{A list, containing any parameters for which you want to -override the default fixed values} +override the default values} } \value{ A transform function diff --git a/man/get_compartment_indices.Rd b/man/get_compartment_indices.Rd new file mode 100644 index 0000000..f716230 --- /dev/null +++ b/man/get_compartment_indices.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/parameters.R +\name{get_compartment_indices} +\alias{get_compartment_indices} +\title{Get compartment indices used throughout the model} +\usage{ +get_compartment_indices() +} +\value{ +A list containing entries: \code{dim} giving the dimensions of each +compartment, currently: group = 18, vax = 4; +\code{group}, a named list of the array indices corresponding to the +first dimension of model compartments: 16 5-year age bands + SW + PBS; +and \code{vax}, a named list of the vaccine strata used in +the second dimension of model compartments: +\enumerate{ +\item historic smallpox; 2. unvaccinated; 3. one-dose; 4. two-dose. +} + +Use this function wherever you need to refer to these standards in your code +to avoid duplication, and make modifying universally easier. +} +\description{ +A function that gets the compartment indices used in the model +} diff --git a/man/parameters_fixed.Rd b/man/parameters_fixed.Rd index 09c9a00..def47d5 100644 --- a/man/parameters_fixed.Rd +++ b/man/parameters_fixed.Rd @@ -4,7 +4,12 @@ \alias{parameters_fixed} \title{Get fixed parameters for use in the model} \usage{ -parameters_fixed(region, initial_infections, overrides = list()) +parameters_fixed( + region, + initial_infections, + use_ve_D = FALSE, + overrides = list() +) } \arguments{ \item{region}{The region for the parameters, must be either \code{"equateur"} or @@ -12,6 +17,9 @@ parameters_fixed(region, initial_infections, overrides = list()) \item{initial_infections}{The initial number of infections} +\item{use_ve_D}{logical, indicating whether model should allow for vaccine +efficacy against death (above and beyond protection against infection)} + \item{overrides}{A list, containing any parameters for which you want to override the default values} } From 4365dc2bc7590c32b9689f864af718c48a05ab3d Mon Sep 17 00:00:00 2001 From: Lilith Whittles Date: Fri, 11 Oct 2024 16:31:55 +0100 Subject: [PATCH 6/7] fix up tests and simplify helper function --- tests/testthat/helper-mpoxseir.R | 14 +- tests/testthat/test-model-targeted-vax.R | 158 ++++++++++------------- 2 files changed, 73 insertions(+), 99 deletions(-) diff --git a/tests/testthat/helper-mpoxseir.R b/tests/testthat/helper-mpoxseir.R index 608eb10..aab214d 100644 --- a/tests/testthat/helper-mpoxseir.R +++ b/tests/testthat/helper-mpoxseir.R @@ -1,18 +1,19 @@ -reference_pars_targeted_vax <- function() { - pars <- parameters_fixed("equateur", initial_infections = 10) +reference_pars_targeted_vax <- function(region = "equateur", uptake = 0) { + pars <- parameters_fixed(region, initial_infections = 10) n_group <- pars$n_group n_vax <- pars$n_vax + idx <- get_compartment_indices() pars$beta_h <- pars$beta_s <- 0.2 / 12.11 pars$beta_z <- rep(0.4 / 12.11, n_group) - vaccination_campaign_length <- 10 - daily_doses <- matrix(0, nrow = vaccination_campaign_length, ncol = n_vax) + pars$vaccination_campaign_length <- 10 + daily_doses <- matrix(0, nrow = pars$vaccination_campaign_length, ncol = n_vax) # nothing happens for j=1 and 2 doses are the maximum - idx <- get_compartment_indices() + daily_doses[, idx$vax$unvaccinated] <- daily_doses[, idx$vax$one_dose] <- 1000 - daily_doses[vaccination_campaign_length, ] <- 0 + daily_doses[pars$vaccination_campaign_length, ] <- 0 pars$daily_doses <- daily_doses pars$N_prioritisation_steps <- 3 @@ -21,6 +22,7 @@ reference_pars_targeted_vax <- function() { c(rep(1,n_group))) # all pars$vaccination_coverage_target_1st_dose_prop <- 0.8 pars$vaccination_coverage_target_2nd_dose_prop <- 0.5 + pars$vaccine_uptake[] <- uptake pars } diff --git a/tests/testthat/test-model-targeted-vax.R b/tests/testthat/test-model-targeted-vax.R index 863cd75..55f1b4d 100644 --- a/tests/testthat/test-model-targeted-vax.R +++ b/tests/testthat/test-model-targeted-vax.R @@ -10,7 +10,7 @@ test_that("run is equal to reference", { expect_true(any(res["cases_inc", , ] > 0)) expect_true(any(res["deaths_inc", , ] > 0)) - expect_equal(sum(res["N_tot", , ] - sum(pars$N)), 0) + expect_equal(sum(res["N_tot", , ] - sum(pars$N0)), 0) }) @@ -26,7 +26,7 @@ test_that("when beta_h = beta_z = beta_s = 0 there are no new infections", { rownames(res) <- names(unlist(m$info()$index)) expect_true(all(res["cases_inc", , ] == 0)) - expect_equal(sum(res["N_tot", , ] - sum(pars$N)), 0) + expect_equal(sum(res["N_tot", , ] - sum(pars$N0)), 0) }) test_that("when CFR = 0 nobody dies", { @@ -42,7 +42,7 @@ test_that("when CFR = 0 nobody dies", { expect_true(all(res["deaths_inc", , ] == 0)) expect_true(any(res["R_tot", , ] > 0)) expect_true(all(res["D_tot", , ] == 0)) - expect_equal(sum(res["N_tot", , ] - sum(pars$N)), 0) + expect_equal(sum(res["N_tot", , ] - sum(pars$N0)), 0) }) test_that("when CFR = 1 everybody dies", { @@ -58,52 +58,60 @@ test_that("when CFR = 1 everybody dies", { expect_true(any(res["deaths_inc", , ] > 0)) expect_true(all(res["R_tot", , ] == 0)) expect_true(any(res["D_tot", , ] > 0)) - expect_equal(sum(res["N_tot", , ] - sum(pars$N)), 0) + expect_equal(sum(res["N_tot", , ] - sum(pars$N0)), 0) }) -test_that("when beta_h = 0 and beta_s=0 there are only zoonotic infections", { +test_that("when beta_h = 0 and beta_s = 0 there are only zoonotic infections", { pars <- reference_pars_targeted_vax() pars$beta_h <- 0 pars$beta_s <- 0 - pars$beta_z<- c(rep(0,pars$n_group-1),0.4 / 12.11) # last group only for test purpose + pars$beta_z[-pars$n_group] <- 0 # last group only for test purpose (i.e. SW) + n_init <- sum(pars$Ea0) + pars$Ea0[] <- 0 m <- model_targeted_vax$new(pars, 1, 3, seed = 1) t <- seq(1, 21) res <- m$simulate(t) - rownames(res) <- names(unlist(m$info()$index)) - - ## isolate the compartments which can have infections (e.g. the last group) - comps_of_interest_Ea <- paste0("Ea", - seq(from=pars$n_group, - length.out=pars$n_vax,by=pars$n_group)) - comps_of_interest_Eb <- paste0("Eb",seq(from=pars$n_group, - length.out=pars$n_vax,by=pars$n_group)) - - expect_true(all(res[grep("^Ea",rownames(res)),,][!(rownames(res)[grep("^Ea",rownames(res))] %in% comps_of_interest_Ea)]==0)) - expect_true(all(res[grep("^Eb",rownames(res)),,][!(rownames(res)[grep("^Eb",rownames(res))] %in% comps_of_interest_Eb)]==0)) - - expect_equal(sum(res["N_tot", , ] - sum(pars$N)), 0) - + idx <- m$info()$index + rownames(res) <- names(unlist(idx)) + y <- m$transform_variables(res) + names(y) <- names(idx) + + # check only infections are in SW + expect_equal(y$cases_cumulative_SW, y$cases_cumulative) + expect_equal(sum(y$I[-pars$n_group, , , ]), 0) + expect_equal(sum(y$N_tot + n_init - sum(pars$N0)), 0) }) -test_that("when beta_h = 0 and beta_z=0 infections only from sexual contact", { +test_that("when beta_h = 0 and beta_z = 0 infections only from sexual contact", { pars <- reference_pars_targeted_vax() pars$beta_h <- 0 - pars$beta_z<- rep(0,pars$n_group) + pars$beta_z[] <- 0 + pars$beta_s <- 0.2 + pars$m_sex["SW", "PBS"] <- pars$m_sex["PBS", "SW"] <- 0.5 + # need to seed infections as so low otherwise - pars$Ir0[17:18,2] <- pars$Ir0[17:18,2] + 1 - pars$Id0[17:18,2] <- pars$Id0[17:18,2] + 1 - pars$S0[17:18,2] <- pars$S0[17:18,2] - pars$Id0[17:18,2] - pars$Ir0[17:18,2] + + idx_comp <- get_compartment_indices() + idx_kp <- unlist(idx_comp$group[c("SW", "PBS")]) + idx_unvax <- idx_comp$vax$unvaccinated + + pars$Ir0[idx_kp, idx_unvax] <- pars$Ir0[idx_kp, idx_unvax] + 10 + pars$Id0[idx_kp, idx_unvax] <- pars$Id0[idx_kp, idx_unvax] + 10 + pars$S0[idx_kp, idx_unvax] <- pars$S0[idx_kp, idx_unvax] - + pars$Id0[idx_kp, idx_unvax] - pars$Ir0[idx_kp, idx_unvax] m <- model_targeted_vax$new(pars, 1, 3, seed = 1) - t <- seq(1, 50) # run for longer to ensure infections in PBS + t <- seq(1, 21) res <- m$simulate(t) - rownames(res) <- names(unlist(m$info()$index)) + idx <- m$info()$index + rownames(res) <- names(unlist(idx)) + ## should have cases in CSW and PBS - expect_true(all(res["cases_cumulative_PBS",,max(t)]>0)) + expect_true(all(res["cases_cumulative_PBS", , max(t)]>0)) expect_true(all(res["cases_cumulative_SW",,max(t)]>0)) ## shouldn't have cases in the age groups @@ -112,7 +120,7 @@ test_that("when beta_h = 0 and beta_z=0 infections only from sexual contact", { expect_true(all(res["cases_cumulative_15_plus",,max(t)]==0)) ## make sure population size continues behaving - expect_equal(sum(res["N_tot", , ] - sum(pars$N)), 0) + expect_equal(sum(res["N_tot", , ] - sum(pars$N0)), 0) }) @@ -120,78 +128,43 @@ test_that("when beta_h = 0 and beta_z=0 infections only from sexual contact", { test_that("when n_vaccination>0, vaccinations are given", { - pars <- reference_pars_targeted_vax() + pars <- reference_pars_targeted_vax(uptake = 0.8) + pars$prioritisation_strategy[] <- 1 + m <- model_targeted_vax$new(pars, 1, 3, seed = 1) t <- seq(1, 21) res <- m$simulate(t) - rownames(res) <- names(unlist(m$info()$index)) + idx <- m$info()$index + rownames(res) <- names(unlist(idx)) + y <- m$transform_variables(res) + names(y) <- names(idx) + idx_comp <- get_compartment_indices() + idx_vax <- c(idx_comp$vax$one_dose, idx_comp$vax$two_dose) - comps_of_interest_S <- paste0("S", - seq(from=2*pars$n_group+1, - to=pars$n_group*pars$n_vax,by=1)) - comps_of_interest_Ea <- paste0("Ea", - seq(from=2*pars$n_group+1, - to=pars$n_group*pars$n_vax,by=1)) - comps_of_interest_Eb <- paste0("Eb", - seq(from=2*pars$n_group+1, - to=pars$n_group*pars$n_vax,by=1)) - comps_of_interest_R <- paste0("R", - seq(from=2*pars$n_group+1, - to=pars$n_group*pars$n_vax,by=1)) ## in first time point there should be no one vaccinated - expect_true(all(res[comps_of_interest_S, ,1] == 0)) - expect_true(all(res[comps_of_interest_Ea, ,1] == 0)) - expect_true(all(res[comps_of_interest_Eb, ,1] == 0)) - expect_true(all(res[comps_of_interest_R, ,1] == 0)) - - ## for subsequent time points we should have people in these classes - expect_true(any(res[comps_of_interest_S,,]>0)) - expect_true(any(res[comps_of_interest_Ea,,]>0)) - expect_true(any(res[comps_of_interest_Eb,,]>0)) - expect_true(any(res[comps_of_interest_R,,]>0)) + expect_true(all(y$N[, idx_vax, , 1] == 0)) + # at last time point we should have lots vaccinated + expect_true(all(y$N[, idx_vax, , 21] > 0)) ## vaccines given should be positive expect_true(any(res[grep("vax_given",rownames(res)),,]>0)) ## population check - expect_equal(sum(res["N_tot", , ] - sum(pars$N)), 0) - -}) - - -test_that("when time is after vaccination_campaign_length, no vaccines are given", { - pars <- reference_pars_targeted_vax() - - m <- model_targeted_vax$new(pars, 1, 3, seed = 1) - t <- seq(1, 21) - res <- m$simulate(t) - rownames(res) <- names(unlist(m$info()$index)) - - expect_true(all(res[grep("vax_given",rownames(res)),,seq(pars$vaccination_campaign_length+1,max(t),1)]==0)) - + expect_equal(sum(res["N_tot", , ] - sum(pars$N0)), 0) + + ## when time is after vaccination_campaign_length, no vaccines are given + expect_equal(sum(y$vax_given_S[, , t > pars$vaccination_campaign_length]), 0) + expect_true(all(y$vax_given_S[, , t == pars$vaccination_campaign_length] > 0)) + + ## the vaccines given do not exceed the total set out in the strategy + expect_true(max(res["total_vax", , ]) <= sum(pars$daily_doses)) }) -test_that("the vaccines given do not exceed the total set out in the strategy", { - pars <- reference_pars_targeted_vax() - - m <- model_targeted_vax$new(pars, 1, 3, seed = 1) - t <- seq(1, 21) - res <- m$simulate(t) - rownames(res) <- names(unlist(m$info()$index)) - - # vaccines_planned <- sum(pars$n_vaccination) - # vaccines_given <- res["total_vax",,max(t)] - - expect_true(all((sum(pars$daily_doses) - res["total_vax",,max(t)])>0)) - - -}) - test_that("if prioritisation_step==1, vaccines are only given in groups 1 - 3", { pars <- reference_pars_targeted_vax() @@ -224,8 +197,7 @@ test_that("if prioritisation_step==1, vaccines are only given in groups 1 - 3", }) test_that("if vaccine_uptake = 0.5, half the expected vaccines are given out", { - pars <- reference_pars_targeted_vax() - pars$vaccine_uptake <- pars$vaccine_uptake*0.5 + pars <- reference_pars_targeted_vax(uptake = 0.5) m <- model_targeted_vax$new(pars, 1, 3, seed = 1) t <- seq(1, 21) @@ -233,14 +205,14 @@ test_that("if vaccine_uptake = 0.5, half the expected vaccines are given out", { rownames(res) <- names(unlist(m$info()$index)) # doses given out should be roughly half of those allocated (roughly comes from rounding due to lack of multinomial) - expect_true(all(res["total_vax",,max(t)]<=0.5*sum(pars$daily_doses))) + expect_true(all(res["total_vax", , max(t)]<=0.5*sum(pars$daily_doses))) }) test_that("no one moves in or out of j=1 (previous smallpox vaccine)", { - pars <- reference_pars_targeted_vax() + pars <- reference_pars_targeted_vax(uptake = 1) m <- model_targeted_vax$new(pars, 1, 3, seed = 1) t <- seq(1, 21) @@ -251,13 +223,13 @@ test_that("no one moves in or out of j=1 (previous smallpox vaccine)", { expect_true(all(res[paste0("N",seq(1:pars$n_group)),,]==pars$S0[,1])) ## population check - expect_equal(sum(res["N_tot", , ] - sum(pars$N)), 0) + expect_equal(sum(res["N_tot", , ] - sum(pars$N0)), 0) }) test_that("1st and 2nd doses are given", { - pars <- reference_pars_targeted_vax() + pars <- reference_pars_targeted_vax(uptake = 1) m <- model_targeted_vax$new(pars, 1, 3, seed = 1) t <- seq(1, 21) @@ -274,7 +246,7 @@ test_that("1st and 2nd doses are given", { test_that("1st and 2nd prioritisation steps can be different depending on the targets", { - pars <- reference_pars_targeted_vax() + pars <- reference_pars_targeted_vax(uptake = 1) # low vaccination target coverage pars$vaccination_coverage_target_1st_dose_prop <- 0.1 pars$vaccination_campaign_length <- 15 @@ -298,7 +270,7 @@ test_that("1st and 2nd prioritisation steps can be different depending on the ta test_that("no 2nd doses are given if no 1st doses are given", { - pars <- reference_pars_targeted_vax() + pars <- reference_pars_targeted_vax(uptake = 1) # no 1st doses allocated pars$daily_doses[,2] <- 0 @@ -317,7 +289,7 @@ test_that("no 2nd doses are given if no 1st doses are given", { test_that("Test compiled compare components", { - pars <- reference_pars_targeted_vax() + pars <- reference_pars_targeted_vax(uptake = 1) pars$exp_noise <- Inf nms <- reference_names() From 857e17d4db0938637c2ad6e62a3bda62f3cb496e Mon Sep 17 00:00:00 2001 From: Lilith Whittles Date: Fri, 11 Oct 2024 17:07:20 +0100 Subject: [PATCH 7/7] bump version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index b9b27cf..8b26440 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: mpoxseir Title: Stochastic compartmental model of mpox transmission -Version: 0.1.1 +Version: 0.1.2 Authors@R: c(person("Lilith", "Whittles", role = c("aut", "cre"), email = "l.whittles@imperial.ac.uk"), person("Ruth", "McCabe", role = c("aut")),