diff --git a/R/disease_progression.R b/R/disease_progression.R index 05367252..fac51866 100644 --- a/R/disease_progression.R +++ b/R/disease_progression.R @@ -84,3 +84,37 @@ update_to_asymptomatic_infection <- function( ) } } + +#' @title Sub-patent progression for vivax model +#' @description +#' randomly selects individuals to clear sub-patent infections based on impact of +#' individual immunity to detectability +#' +#' @param state the human state variable +#' @param variables the available human variables +#' @param parameters model parameters +#' @noRd +create_subpatent_progression_process <- function( + state, + variables, + parameters + ){ + function(timestep) { + ## Get recovery rates for each subpatent individual + subpatent_index <- state$get_index_of("U") + rate <- anti_parasite_immunity( + parameters$du_min, parameters$du_max, parameters$au50, parameters$ku, + variables$id$get_values(index = subpatent_index), + variables$idm$get_values(index = subpatent_index)) + ## Get individuals who are going to change + to_move <- subpatent_index$sample(1/rate) + # update individuals and infectivity + update_infection( + state, + to_state = "S", + infectivity = variables$infectivity, + new_infectivity = 0, + to_move + ) + } +} diff --git a/R/human_infection.R b/R/human_infection.R index 8d5ea260..4561cb00 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -53,14 +53,33 @@ simulate_infection <- function( parameters$ud ) } + + if(parameters$parasite == "falciparum"){ + clinical_infections <- calculate_clinical_infections( + variables, + infected_humans, + parameters, + renderer, + timestep + ) + } else if (parameters$parasite == "vivax"){ + patent_infections <- calculate_patent_infections( + variables, + infected_humans, + parameters, + renderer, + timestep + ) + + clinical_infections <- calculate_clinical_infections( + variables, + patent_infections, + parameters, + renderer, + timestep + ) + } - clinical_infections <- calculate_clinical_infections( - variables, - infected_humans, - parameters, - renderer, - timestep - ) update_severe_disease( timestep, @@ -175,15 +194,51 @@ calculate_infections <- function( infected } -#' @title Calculate clinical infections +#' @title Calculate patent infections (vivax only) #' @description -#' Sample clinical infections from all infections +#' Sample patent infections from all infections #' @param variables a list of all of the model variables #' @param infections bitset of infected humans #' @param parameters model parameters #' @param renderer model render #' @param timestep current timestep #' @noRd +calculate_patent_infections <- function( + variables, + infections, + parameters, + renderer, + timestep +) { + id <- variables$id$get_values(infections) + idm <- variables$idm$get_values(infections) + philm <- anti_parasite_immunity( + min = parameters$du_min, max = parameters$du_max, a50 = parameters$au50, + k = parameters$ku, id = id, idm = idm) + patent_infections <- bitset_at(infections, bernoulli_multi_p(philm)) + incidence_renderer( + variables$birth, + renderer, + patent_infections, + infections, + philm, + 'inc_patent_', + parameters$patent_incidence_rendering_min_ages, + parameters$patent_incidence_rendering_max_ages, + timestep + ) + patent_infections +} + +#' @title Calculate clinical infections +#' @description +#' Sample clinical infections from all infections or patent infections (vivax) +#' @param variables a list of all of the model variables +#' @param infections bitset of infected/patent humans +#' @param parameters model parameters +#' @param renderer model render +#' @param timestep current timestep +#' @noRd calculate_clinical_infections <- function( variables, infections, @@ -458,3 +513,9 @@ blood_immunity <- function(ib, parameters) { (1 + (ib / parameters$ib0) ** parameters$kb) ) } + +# Implemented from White et al., 2018 - Supplementary Information +anti_parasite_immunity <- function(min, max, a50, k, id, idm){ + min + (max - min) / ( + 1 + ((id + idm) / a50) ** k) +} diff --git a/R/model.R b/R/model.R index 48af9af9..3c83172b 100644 --- a/R/model.R +++ b/R/model.R @@ -29,22 +29,28 @@ #' * n: number of humans between an inclusive age range at this timestep. This #' defaults to n_730_3650. Other age ranges can be set with #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters. -#' * n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This -#' defaults to n_detect_730_3650. Other age ranges can be set with +#' * n_detect_pcr: number of humans with an infection detectable by pcr between an inclusive age range at this timestep. This +#' defaults to n_detect_pcr_730_3650. Other age ranges can be set with #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters. -#' * p_detect: the sum of probabilities of detection by microscopy between an +#' * n_detect_lm: number of humans with an infection detectable by light microscopy between an inclusive age range at this timestep. This +#' defaults to n_detect_lm_730_3650. Other age ranges can be set with +#' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters. +#' * p_detect_lm: the sum of probabilities of detection by light microscopy between an #' inclusive age range at this timestep. This -#' defaults to p_detect_730_3650. Other age ranges can be set with +#' defaults to p_detect_lm_730_3650. Other age ranges can be set with #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters. -#' * n_severe: number of humans with a severe infection between an inclusive -#' age range at this timestep. Age ranges can be set with -#' severe_prevalence_rendering_min_ages and severe_prevalence_rendering_max_ages parameters. #' * n_inc: number of new infections for humans between an inclusive age range at this timestep. #' incidence columns can be set with #' incidence_rendering_min_ages and incidence_rendering_max_ages parameters. #' * p_inc: sum of probabilities of infection for humans between an inclusive age range at this timestep. #' incidence columns can be set with #' incidence_rendering_min_ages and incidence_rendering_max_ages parameters. +#' * n_inc_patent: number of new patent infections for humans between an inclusive age range at this timestep. +#' patent incidence columns can be set with +#' patent_incidence_rendering_min_ages and patent_incidence_rendering_max_ages parameters. +#' * p_inc_patent: sub of probabilities of patent infection for humans between an inclusive age range at this timestep. +#' patent incidence columns can be set with +#' patent_incidence_rendering_min_ages and patent_incidence_rendering_max_ages parameters. #' * n_inc_clinical: number of new clinical infections for humans between an inclusive age range at this timestep. #' clinical incidence columns can be set with #' clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters. diff --git a/R/mortality_processes.R b/R/mortality_processes.R index 65b090cb..3d406fd6 100644 --- a/R/mortality_processes.R +++ b/R/mortality_processes.R @@ -69,6 +69,11 @@ sample_maternal_immunity <- function(variables, target, timestep, parameters) { birth_ivm <- variables$ica$get_values(mothers) * parameters$pvm variables$icm$queue_update(birth_icm, target_group) variables$ivm$queue_update(birth_ivm, target_group) + + if(parameters$parasite == "vivax"){ + birth_idm <- variables$id$get_values(mothers) * parameters$pcm + variables$idm$queue_update(birth_idm, target_group) + } } } } diff --git a/R/parameters.R b/R/parameters.R index e611464c..9abc748b 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -5,14 +5,14 @@ #' are for Plasmodium falciparum, are explained in #' "The US President's Malaria Initiative, Plasmodium falciparum transmission #' and mortality: A modelling study." -#' +#' #' Plasmodium vivax specific parameters are explained in #' "Mathematical modelling of the impact of expanding levels of malaria control #' interventions on Plasmodium vivax." by White, Michael T., et al. #' #' @param overrides a named list of parameter values to use instead of defaults #' @param parasite Plasmodium parasite species ("falciparum" or "vivax"); default = "falciparum" -#' +#' #' The parameters are defined below. #' #' human fixed state transitions: @@ -23,7 +23,7 @@ #' * du - the delay for humans to move from state U to S (P.f only); default = 110 #' #' duration of subpatent infection: du (P.v only): -#' +#' #' * du_max - Maximum duration of subpatent infection: default = 70 #' * du_min - Minimum duration of subpatent infection: default = 10 #' * ku - Shape parameter: default = 4.602 @@ -37,7 +37,7 @@ #' * init_ivm - the immunity from severe disease at birth (P.f only); default = 0 #' * init_ib - the initial pre-erythrocitic immunity (P.f only); default = 0 #' * init_id - the initial acquired immunity to detectability; default = 0 -#' +#' #' immunity boost grace periods: #' #' * ub - period in which pre-erythrocytic immunity cannot be boosted (P.f only); default = 7.2 @@ -46,10 +46,10 @@ #' * uv - period in which severe immunity cannot be boosted (P.f only); default = 11.4321 #' #' maternal immunity parameters: -#' +#' #' * pcm - new-born clinical immunity relative to mother's; default = 0.774368 #' * pvm - new-born severe immunity relative to mother's (P.f only); default = 0.195768 -#' +#' #' immunity decay rates: #' #' * rm - decay rate for maternal immunity to clinical disease; default = 67.6952 @@ -67,7 +67,7 @@ #' * kb - shape parameter (P.f only); default = 2.16 #' * b - probability of pre-erythrocytic infection (P.v only): default = 0.5 #' -#' probability of asymptomatic detection/detectbility immunity (P.f only): +#' probability of asymptomatic detection/detectability immunity (P.f only): #' #' * fd0 - time-scale at which immunity changes with age; default = 0.007055 #' * ad - scale parameter relating age to immunity; default = 7993.5 @@ -82,7 +82,7 @@ #' * phi1lm - maximum reduction due to immunity; default = 0.00482170890334156 #' * ic0lm - scale parameter; default = 27.52 #' * kclm - shape parameter; default = 2.403 -#' +#' #' probability of clinical infection/clinical immunity: #' #' * phi0 - maximum probability due to no immunity; default = 0.792 @@ -113,16 +113,16 @@ #' * de - Duration of the human latent period of infection; default = 12 #' * delay_gam - Lag from parasites to infectious gametocytes; default = 12.5 #' * dem - Extrinsic incubation period in mosquito population model; default = 10 -#' +#' #' hypnozoite parameters (P.v only): -#' +#' #' * f - relapse rate; default = 0.024 #' * gammal - clearance rate; default = 0.0026 -#' +#' #' parasite parameter -#' +#' #' * parasite - parasite species (falciparum or vivax); default = "falciparum" -#' +#' #' human population parameters: #' * human_population - the initial number of humans to model; default = 100 #' * human_population_timesteps - the timesteps at which the population should @@ -130,7 +130,7 @@ #' * custom_demography - population demography given; default = FALSE, #' * average_age - the average age of humans (in timesteps), this is only used #' if custom_demography is FALSE; default = 7663 -#' +#' #' initial state proportions: #' #' * s_proportion - the proportion of `human_population` that begin as susceptible; default = 0.420433246 @@ -143,25 +143,25 @@ #' * t_proportion - the proportion of `human_population` that begin treated; default = 0 #' #' rendering: -#' All values are in timesteps and all ranges are inclusive -#' please set rendered age groups using the convenience function +#' All values are in timesteps and all ranges are inclusive. +#' Please set rendered age groups using the convenience function #' `set_demography` #' -#' * age_group_rendering_min_ages - the minimum ages for population size outputs; default = numeric(0) -#' * age_group_rendering_max_ages - the corresponding max ages; default = numeric(0) +#' * age_group_rendering_min_ages - the minimum ages for population size outputs; default = turned off +#' * age_group_rendering_max_ages - the corresponding max ages; default = turned off #' * prevalence_rendering_min_ages - the minimum ages for clinical prevalence -#' outputs; default = 730 +#' outputs (pcr and lm detectable infections); default = 730 #' * prevalence_rendering_max_ages - the corresponding max ages; default = 3650 #' * incidence_rendering_min_ages - the minimum ages for incidence -#' outputs (includes asymptomatic microscopy +); default = turned off +#' outputs (pf: D+Tr+A, pv: D+Tr+A+U); default = turned off #' * incidence_rendering_max_ages - the corresponding max ages; default = turned off -#' * clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = 0 -#' * clinical_incidence_rendering_max_ages - the corresponding max ages; default = 1825 +#' * patent_incidence_rendering_min_ages - the minimum ages for patent incidence outputs (lm detectable), (P.v only); default = turned off +#' * patent_incidence_rendering_max_ages - the corresponding max ages (P.v only); default = turned off +#' * clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = turned off +#' * clinical_incidence_rendering_max_ages - the corresponding max ages; default = turned off #' * severe_incidence_rendering_min_ages - the minimum ages for severe incidence #' outputs; default = turned off #' * severe_incidence_rendering_max_ages - the corresponding max ages; default = turned off -#' * severe_prevalence_rendering_min_ages - the minimum ages for severe prevalance outputs; default = numeric(0), -#' * severe_prevalence_rendering_max_ages - the corresponding max ages; default = numeric(0), #' #' mosquito life stage transitions: #' @@ -192,11 +192,11 @@ #' flexible carrying capacity parameters: #' please set carrying capacity using the convenience function #' `set_carrying_capacity` -#' +#' #' * carrying_capacity; default = FALSE #' * carrying_capacity_timesteps; default = NULL #' * carrying_capacity_values; default = NULL -#' +#' #' larval mortality rates: #' #' * me - early stage larval mortality rate; default = 0.0338 @@ -241,7 +241,7 @@ #' clinical treatment coverage changes (these values refer to the index in drug_* parameters); default = NULL #' * clinical_treatment_coverage - a list of vectors giving coverage values for each drug; default = NULL #' -#' PEV parameters: +#' PEV parameters: #' please set vaccine strategies with the convenience functions #' `set_pev_epi` and `set_mass_pev` #' @@ -278,7 +278,7 @@ #' #' #' miscellaneous: -#' +#' #' * mosquito_limit - the maximum number of mosquitoes to allow for in the #' simulation; default = 1.00E+05 #' * individual_mosquitoes - boolean whether adult mosquitoes are modeled @@ -292,7 +292,7 @@ #' #' @export get_parameters <- function(overrides = list(), parasite = "falciparum") { - + parameters <- c( ## Parasite-specific parameters set in parasite_parameters.csv # human fixed state transitions @@ -302,7 +302,7 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") { # maternal immunity parameters # immunity decay rates # probability of pre-erythrocytic infection/blood immunity - # probability of asymptomatic detection/detectbility immunity + # probability of asymptomatic detection/detectability immunity # probability of asymptomatic infection/antiparasite immunity # probability of clinical infection/clinical immunity # probability of severe infection/severe disease immunity @@ -310,12 +310,12 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") { # parasite incubation periods # hypnozoite parameters parasite_parameters[[parasite]], - + ## Add remaining parameters list( # parasite parameter parasite = parasite, - + # human population parameters human_population = 100, human_population_timesteps = 0, @@ -334,13 +334,15 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") { prevalence_rendering_max_ages = 10 * 365, incidence_rendering_min_ages = numeric(0), incidence_rendering_max_ages = numeric(0), + patent_incidence_rendering_min_ages = numeric(0), + patent_incidence_rendering_max_ages = numeric(0), clinical_incidence_rendering_min_ages = numeric(0), - clinical_incidence_rendering_max_ages = 5 * 365, + clinical_incidence_rendering_max_ages = numeric(0), severe_incidence_rendering_min_ages = numeric(0), severe_incidence_rendering_max_ages = numeric(0), severe_prevalence_rendering_min_ages = numeric(0), severe_prevalence_rendering_max_ages = numeric(0), - + # mosquito life stage transitions del = 6.64, dl = 3.72, @@ -430,7 +432,7 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") { tbv_timesteps = NULL, tbv_coverages = NULL, tbv_ages = NULL, - + # misc mosquito_limit = 100 * 1000, individual_mosquitoes = FALSE, @@ -440,19 +442,19 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") { ode_max_steps = 1e6, progress_bar = FALSE )) - + # Override parameters with any client specified ones if (!is.list(overrides)) { stop('overrides must be a list') } - + for (name in names(overrides)) { if (!(name %in% names(parameters))) { stop(paste('unknown parameter', name, sep=' ')) } parameters[[name]] <- overrides[[name]] } - + props <- c( parameters$s_proportion, parameters$d_proportion, @@ -460,11 +462,11 @@ get_parameters <- function(overrides = list(), parasite = "falciparum") { parameters$u_proportion, parameters$t_proportion ) - + if (!approx_sum(props, 1)) { stop("Starting proportions do not sum to 1") } - + parameters } @@ -539,7 +541,7 @@ parameterise_total_M <- function(parameters, total_M) { } #' Use parameter draw from the join posterior -#' +#' #' Overrides default (median) model parameters with a single draw from the fitted #' joint posterior. Must be called prior to set_equilibrium. #' diff --git a/R/processes.R b/R/processes.R index 2b431122..da62096f 100644 --- a/R/processes.R +++ b/R/processes.R @@ -46,7 +46,14 @@ create_processes <- function( create_exponential_decay_process(variables$iva, parameters$rva), create_exponential_decay_process(variables$id, parameters$rid) ) - + + if(parameters$parasite == "vivax"){ + processes <- c( + processes, + # Maternal immunity to detectability + create_exponential_decay_process(variables$idm, parameters$rm)) + } + if (parameters$individual_mosquitoes) { processes <- c( processes, @@ -59,7 +66,7 @@ create_processes <- function( ) ) } - + # ============================== # Biting and mortality processes # ============================== @@ -96,7 +103,7 @@ create_processes <- function( parameters$da, variables$infectivity, parameters$ca - ) + ) }, create_progression_process( variables$state, @@ -106,14 +113,21 @@ create_processes <- function( variables$infectivity, parameters$cu ), - create_progression_process( - variables$state, - 'U', - 'S', - parameters$du, - variables$infectivity, - 0 - ), + if(parameters$parasite == "falciparum"){ + create_progression_process( + variables$state, + 'U', + 'S', + parameters$du, + variables$infectivity, + 0) + } else if (parameters$parasite == "vivax"){ + create_subpatent_progression_process( + variables$state, + variables, + parameters + ) + }, create_progression_process( variables$state, 'Tr', @@ -123,7 +137,7 @@ create_processes <- function( 0 ) ) - + # =============== # ODE integration # =============== @@ -131,7 +145,7 @@ create_processes <- function( processes, create_solver_stepping_process(solvers, parameters) ) - + # ========= # RTS,S EPI # ========= @@ -148,7 +162,7 @@ create_processes <- function( ) ) } - + # ========= # PMC # ========= @@ -167,21 +181,27 @@ create_processes <- function( ) ) } - + # ========= # Rendering # ========= + + imm_var_names <- c('ica','icm','ib','id','iva','ivm') + if(parameters$parasite == "vivax"){ + imm_var_names <- c(imm_var_names, 'idm') + } + processes <- c( processes, individual::categorical_count_renderer_process( renderer, variables$state, - c('S', 'A', 'D', 'U', 'Tr') + c('S', 'D', 'A', 'U', 'Tr') ), 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( variables$state, @@ -197,7 +217,7 @@ create_processes <- function( ), create_compartmental_rendering_process(renderer, solvers, parameters) ) - + if (parameters$individual_mosquitoes) { processes <- c( processes, @@ -219,11 +239,11 @@ create_processes <- function( ) ) } - + # ====================== # Intervention processes # ====================== - + if (parameters$bednets) { processes <- c( processes, @@ -236,14 +256,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 # ====================== @@ -253,7 +273,7 @@ create_processes <- function( create_progress_process(timesteps) ) } - + processes } @@ -277,7 +297,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) @@ -294,7 +314,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 index 3fcf4882..70fa0937 100644 --- a/R/render.R +++ b/R/render.R @@ -3,65 +3,91 @@ 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 -#' +#' +#' @description renders prevalence numerators and denominators for individuals +#' detected by lm microscopy and pcr, where those infected asymptomatically by +#' P. falciparum have reduced probability of infection due to detectability +#' immunity (reported as an integer sample (n_detect_lm_) and summing over +#' detection probabilities: p_detect_lm) +#' #' @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( - state, - birth, - immunity, - parameters, - renderer - ) { + state, + birth, + immunity, + parameters, + renderer +) { function(timestep) { + asymptomatic <- state$get_index_of('A') - prob <- probability_of_detection( - get_age(birth$get_values(asymptomatic), timestep), - immunity$get_values(asymptomatic), - parameters - ) - asymptomatic_detected <- bitset_at(asymptomatic, bernoulli_multi_p(prob)) + + if(parameters$parasite == "falciparum"){ + prob <- probability_of_detection( + get_age(birth$get_values(asymptomatic), timestep), + immunity$get_values(asymptomatic), + parameters + ) + asymptomatic_detected <- bitset_at(asymptomatic, bernoulli_multi_p(prob)) + + } else { + asymptomatic_detected <- asymptomatic + } clinically_detected <- state$get_index_of(c('Tr', 'D')) - detected <- clinically_detected$copy()$or(asymptomatic_detected) + lm_detected <- clinically_detected$or(asymptomatic_detected) + pcr_detected <- state$get_index_of(c('Tr', 'D', 'A', 'U')) for (i in seq_along(parameters$prevalence_rendering_min_ages)) { lower <- parameters$prevalence_rendering_min_ages[[i]] upper <- parameters$prevalence_rendering_max_ages[[i]] in_age <- in_age_range(birth, timestep, lower, upper) + + # render age renderer$render( paste0('n_', lower, '_', upper), in_age$size(), timestep - ) + ) + + # render pcr detection renderer$render( - paste0('n_detect_', lower, '_', upper), - in_age$copy()$and(detected)$size(), + paste0('n_detect_pcr_', lower, '_', upper), + pcr_detected$and(in_age)$size(), timestep ) + + # render lm detection renderer$render( - paste0('p_detect_', lower, '_', upper), - in_age$copy()$and(clinically_detected)$size() + sum( - prob[bitset_index(asymptomatic, in_age)] - ), + paste0('n_detect_lm_', lower, '_', upper), + lm_detected$and(in_age)$size(), timestep ) + + if(parameters$parasite == "falciparum"){ + # render lm detection (falciparum): summed probability + renderer$render( + paste0('p_detect_lm_', lower, '_', upper), + clinically_detected$and(in_age)$size() + sum( + prob[bitset_index(asymptomatic, in_age)] + ), + timestep + ) + } } } } #' @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,19 +97,19 @@ create_prevelance_renderer <- function( #' @param lowers age bounds #' @param uppers age bounds #' @param timestep current target -#' +#' #' @noRd incidence_renderer <- function( - birth, - renderer, - target, - source_pop, - prob, - prefix, - lowers, - uppers, - timestep - ) { + birth, + renderer, + target, + source_pop, + prob, + prefix, + lowers, + uppers, + timestep +) { for (i in seq_along(lowers)) { lower <- lowers[[i]] upper <- uppers[[i]] @@ -104,9 +130,9 @@ incidence_renderer <- function( } create_variable_mean_renderer_process <- function( - renderer, - names, - variables + renderer, + names, + variables ) { function(timestep) { for (i in seq_along(variables)) { @@ -120,12 +146,12 @@ create_variable_mean_renderer_process <- function( } create_vector_count_renderer_individual <- function( - mosquito_state, - species, - state, - renderer, - parameters - ) { + mosquito_state, + species, + state, + renderer, + parameters +) { function(timestep) { adult <- mosquito_state$get_index_of('NonExistent')$not(TRUE) for (i in seq_along(parameters$species)) { @@ -168,7 +194,7 @@ create_age_group_renderer <- function( paste0('n_age_', lower, '_', upper), in_age$size(), timestep - ) + ) } } } diff --git a/R/variables.R b/R/variables.R index 860b3396..bcb74669 100644 --- a/R/variables.R +++ b/R/variables.R @@ -96,6 +96,19 @@ create_variables <- function(parameters) { last_boosted_id <- individual::DoubleVariable$new(rep(-1, size)) # Maternal immunity + if(parameters$parasite == "vivax"){ + idm <- individual::DoubleVariable$new( + initial_immunity( + parameters$init_idm, + initial_age, + groups, + eq, + parameters, + 'IDM' + ) + ) + } + icm <- individual::DoubleVariable$new( initial_immunity( parameters$init_icm, @@ -224,6 +237,10 @@ create_variables <- function(parameters) { net_time = net_time, spray_time = spray_time ) + + if(parameters$parasite == "vivax"){ + variables <- c(variables, idm = idm) + } # Add variables for individual mosquitoes if (parameters$individual_mosquitoes) { diff --git a/data-raw/parasite_parameters.csv b/data-raw/parasite_parameters.csv index ab054dc5..548b2361 100644 --- a/data-raw/parasite_parameters.csv +++ b/data-raw/parasite_parameters.csv @@ -52,7 +52,6 @@ falciparum,dem,10, vivax,dd,5, vivax,dt,1, vivax,da,10, -vivax,du,110,TO BE REMOVED: NOT VIVAX: Job 3: subpatent duration vivax,du_max,70, vivax,du_min,10, vivax,ku,4.602, diff --git a/data/parasite_parameters.rda b/data/parasite_parameters.rda index 29dde8eb..059a87fb 100644 Binary files a/data/parasite_parameters.rda and b/data/parasite_parameters.rda differ diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd index 92ae4b3a..108dceda 100644 --- a/man/get_parameters.Rd +++ b/man/get_parameters.Rd @@ -72,7 +72,7 @@ probability of pre-erythrocytic infection/blood immunity: \item b - probability of pre-erythrocytic infection (P.v only): default = 0.5 } -probability of asymptomatic detection/detectbility immunity (P.f only): +probability of asymptomatic detection/detectability immunity (P.f only): \itemize{ \item fd0 - time-scale at which immunity changes with age; default = 0.007055 \item ad - scale parameter relating age to immunity; default = 7993.5 @@ -159,25 +159,25 @@ subpatents; default = 0.133028023 } rendering: -All values are in timesteps and all ranges are inclusive -please set rendered age groups using the convenience function +All values are in timesteps and all ranges are inclusive. +Please set rendered age groups using the convenience function \code{set_demography} \itemize{ -\item age_group_rendering_min_ages - the minimum ages for population size outputs; default = numeric(0) -\item age_group_rendering_max_ages - the corresponding max ages; default = numeric(0) +\item age_group_rendering_min_ages - the minimum ages for population size outputs; default = turned off +\item age_group_rendering_max_ages - the corresponding max ages; default = turned off \item prevalence_rendering_min_ages - the minimum ages for clinical prevalence -outputs; default = 730 +outputs (pcr and lm detectable infections); default = 730 \item prevalence_rendering_max_ages - the corresponding max ages; default = 3650 \item incidence_rendering_min_ages - the minimum ages for incidence -outputs (includes asymptomatic microscopy +); default = turned off +outputs (pf: D+Tr+A, pv: D+Tr+A+U); default = turned off \item incidence_rendering_max_ages - the corresponding max ages; default = turned off -\item clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = 0 -\item clinical_incidence_rendering_max_ages - the corresponding max ages; default = 1825 +\item patent_incidence_rendering_min_ages - the minimum ages for patent incidence outputs (lm detectable), (P.v only); default = turned off +\item patent_incidence_rendering_max_ages - the corresponding max ages (P.v only); default = turned off +\item clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = turned off +\item clinical_incidence_rendering_max_ages - the corresponding max ages; default = turned off \item severe_incidence_rendering_min_ages - the minimum ages for severe incidence outputs; default = turned off \item severe_incidence_rendering_max_ages - the corresponding max ages; default = turned off -\item severe_prevalence_rendering_min_ages - the minimum ages for severe prevalance outputs; default = numeric(0), -\item severe_prevalence_rendering_max_ages - the corresponding max ages; default = numeric(0), } mosquito life stage transitions: @@ -304,8 +304,6 @@ please set TBV parameters with the convenience functions } miscellaneous: -please set parasite species using the convenience function -\code{set_parasite} \itemize{ \item mosquito_limit - the maximum number of mosquitoes to allow for in the simulation; default = 1.00E+05 diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd index 64cac9c4..2b8319bd 100644 --- a/man/run_simulation.Rd +++ b/man/run_simulation.Rd @@ -45,22 +45,28 @@ the whole population) \item n: number of humans between an inclusive age range at this timestep. This defaults to n_730_3650. Other age ranges can be set with prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters. -\item n_detect: number of humans with an infection detectable by microscopy between an inclusive age range at this timestep. This -defaults to n_detect_730_3650. Other age ranges can be set with +\item n_detect_pcr: number of humans with an infection detectable by pcr between an inclusive age range at this timestep. This +defaults to n_detect_pcr_730_3650. Other age ranges can be set with prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters. -\item p_detect: the sum of probabilities of detection by microscopy between an +\item n_detect_lm: number of humans with an infection detectable by light microscopy between an inclusive age range at this timestep. This +defaults to n_detect_lm_730_3650. Other age ranges can be set with +prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters. +\item p_detect_lm: the sum of probabilities of detection by light microscopy between an inclusive age range at this timestep. This -defaults to p_detect_730_3650. Other age ranges can be set with +defaults to p_detect_lm_730_3650. Other age ranges can be set with prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters. -\item n_severe: number of humans with a severe infection between an inclusive -age range at this timestep. Age ranges can be set with -severe_prevalence_rendering_min_ages and severe_prevalence_rendering_max_ages parameters. \item n_inc: number of new infections for humans between an inclusive age range at this timestep. incidence columns can be set with incidence_rendering_min_ages and incidence_rendering_max_ages parameters. \item p_inc: sum of probabilities of infection for humans between an inclusive age range at this timestep. incidence columns can be set with incidence_rendering_min_ages and incidence_rendering_max_ages parameters. +\item n_inc_patent: number of new patent infections for humans between an inclusive age range at this timestep. +patent incidence columns can be set with +patent_incidence_rendering_min_ages and patent_incidence_rendering_max_ages parameters. +\item p_inc_patent: sub of probabilities of patent infection for humans between an inclusive age range at this timestep. +patent incidence columns can be set with +patent_incidence_rendering_min_ages and patent_incidence_rendering_max_ages parameters. \item n_inc_clinical: number of new clinical infections for humans between an inclusive age range at this timestep. clinical incidence columns can be set with clinical_incidence_rendering_min_ages and clinical_incidence_rendering_max_ages parameters. diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index a8e6033e..b8cf2118 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -220,9 +220,6 @@ test_that('calculate_clinical_infections correctly samples clinically infected', ) immunity_mock <- mockery::mock(c(.2, .3, .4)) - boost_mock <- mockery::mock() - mockery::stub(calculate_clinical_infections, 'boost_immunity', boost_mock) - mockery::stub(calculate_clinical_infections, 'clinical_immunity', immunity_mock) bernoulli_mock <- mockery::mock(c(1, 3)) mockery::stub(calculate_clinical_infections, 'bernoulli_multi_p', bernoulli_mock) diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R index 920f09bb..2a76f7e9 100644 --- a/tests/testthat/test-render.R +++ b/tests/testthat/test-render.R @@ -34,15 +34,23 @@ test_that('that default rendering works', { mockery::expect_args( renderer$render_mock(), 2, - 'n_detect_730_3650', + 'n_detect_pcr_730_3650', + 3, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 3, + 'n_detect_lm_730_3650', 2, timestep ) mockery::expect_args( renderer$render_mock(), - 3, - 'p_detect_730_3650', + 4, + 'p_detect_lm_730_3650', 1.5, timestep ) diff --git a/tests/testthat/test-vivax.R b/tests/testthat/test-vivax.R index 729bcc93..f08594e6 100644 --- a/tests/testthat/test-vivax.R +++ b/tests/testthat/test-vivax.R @@ -18,9 +18,178 @@ test_that('Test parasite = vivax sets parameters$parasite to vivax', { test_that('Test difference between falciparum and vivax parameter lists', { falciparum_parameters <- get_parameters(parasite = "falciparum") vivax_parameters <- get_parameters(parasite = "vivax") - + expect_true(all(names(falciparum_parameters)[!names(falciparum_parameters) %in% names(vivax_parameters)] %in% c("du","rvm","rva","rb","b0","b1","ib0","kb","theta0","theta1","kv","fv0","av","gammav","iv0","fd0","ad","gammad","d1","id0","kd","ub","uv","gamma1","pvm","init_ivm","init_ib","init_iva"))) expect_true(all(names(vivax_parameters[!names(vivax_parameters) %in% names(falciparum_parameters)]) %in% c("du_max","du_min","ku","au50","b","phi0lm","phi1lm","ic0lm","kclm","d0","ca","init_idm","f","gammal"))) }) + +## Test subpatent progression functions +test_that('Test subpatent duration function works ', { + vivax_parameters <- get_parameters(parasite = "vivax", + overrides = list(s_proportion = 0, + d_proportion = 0, + a_proportion = 0, + u_proportion = 1, + t_proportion = 0)) + variables <- create_variables(vivax_parameters) + index <- variables$state$get_index_of("U") + du_min <- vivax_parameters$du_min + du_max <- vivax_parameters$du_max + au50 <- vivax_parameters$au50 + ku <- vivax_parameters$ku + + expect_equal(object = anti_parasite_immunity( + du_min, du_max, au50, ku, + variables$id$get_values(index), + variables$idm$get_values(index)), + expected = rep(du_max,100)) + + ## Change initial values of ID, and IDM, check they are the same + variables$id <- individual::DoubleVariable$new(1:100) + ID_durations <- anti_parasite_immunity( + du_min, du_max, au50, ku, + variables$id$get_values(index), + variables$idm$get_values(index)) + + variables$id <- individual::DoubleVariable$new(rep(0,100)) + variables$idm <- individual::DoubleVariable$new(1:100) + IDM_durations <- anti_parasite_immunity( + du_min, du_max, au50, ku, + variables$id$get_values(index), + variables$idm$get_values(index)) + + expect_equal(object = ID_durations, expected = IDM_durations) + + ## Check convergence to min_du at high immunity + variables$id <- individual::DoubleVariable$new(rep(1E7,100)) + variables$idm <- individual::DoubleVariable$new(rep(0,100)) + expect_equal(object = anti_parasite_immunity( + du_min, du_max, au50, ku, + variables$id$get_values(index), + variables$idm$get_values(index)), + expected = rep(du_min,100), + tolerance = 1E-2) +}) + +test_that('Test default vivax rendering works', { + + timestep <- 0 + year <- 365 + birth <- individual::IntegerVariable$new( + -c(2, 5, 10, 11) * year + ) + vivax_parameters <- get_parameters( + parasite = "vivax") + + renderer <- mock_render(1) + incidence_renderer( + birth, + renderer, + individual::Bitset$new(4)$insert(c(1, 2, 4)), + individual::Bitset$new(4)$insert(seq(4)), + c(.1, .2, .3, .4), + 'inc_patent_', + c(0, 2) * year, + c(5, 10) * year, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 1, + 'n_0_1825', + 2, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 2, + 'n_inc_patent_0_1825', + 2, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 3, + 'p_inc_patent_0_1825', + 0.3, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 4, + 'n_730_3650', + 3, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 5, + 'n_inc_patent_730_3650', + 2, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 6, + 'p_inc_patent_730_3650', + .6, + timestep + ) +}) + +test_that('that vivax patent prevalence rendering works', { + + timestep <- 1 + state <- individual::CategoricalVariable$new( + c('U', 'A', 'D', 'S', 'Tr'), + c('U', 'A', 'D', 'D', 'D', 'S') + ) + birth <- individual::IntegerVariable$new( + -c(3, 4, 5, 1, 11, 6) * 365 + ) + immunity <- individual::DoubleVariable$new(rep(1, 6)) + vivax_parameters <- get_parameters(parasite = "vivax") + renderer <- mock_render(1) + + process <- create_prevelance_renderer( + state, + birth, + immunity, + vivax_parameters, + renderer + ) + + process(timestep) + + mockery::expect_args( + renderer$render_mock(), + 1, + 'n_730_3650', + 4, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 2, + 'n_detect_pcr_730_3650', + 3, + timestep + ) + + mockery::expect_args( + renderer$render_mock(), + 3, + 'n_detect_lm_730_3650', + 2, + timestep + ) +})