diff --git a/NAMESPACE b/NAMESPACE index 4fd1b1da..4ba1e15d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,6 +18,7 @@ export(rtss_profile) export(run_metapop_simulation) export(run_simulation) export(run_simulation_with_repetitions) +export(set_antimalarial_resistance) export(set_bednets) export(set_carrying_capacity) export(set_clinical_treatment) diff --git a/R/antimalarial_resistance.R b/R/antimalarial_resistance.R new file mode 100644 index 00000000..532c042b --- /dev/null +++ b/R/antimalarial_resistance.R @@ -0,0 +1,139 @@ +#' @title Parameterise antimalarial resistance +#' @description +#' Parameterise antimalarial resistance +#' +#' @param parameters the model parameters +#' @param drug the index of the drug which resistance is being set, as set by the set_drugs() function, in the parameter list +#' @param timesteps vector of time steps for each update to resistance proportion and resistance outcome probability +#' @param artemisinin_resistance vector of updates to the proportions of infections that are artemisinin resistant at time t +#' @param partner_drug_resistance vector of updates to the proportions of infections that are partner-drug resistant at time t +#' @param slow_parasite_clearance_prob vector of updates to the proportion of artemisinin-resistant infections that result in early treatment failure +#' @param early_treatment_failure_prob vector of updates to the proportion of artemisinin-resistant infections that result in slow parasite clearance +#' @param late_clinical_failure_prob vector of updates to the proportion of partner-drug-resistant infections that result in late clinical failure +#' @param late_parasitological_prob vector of updates to the proportion of partner-drug-resistant infections that result in late parasitological failure +#' @param reinfection_prob vector of updates to the proportion of partner-drug-resistant infections that result in reinfection during prophylaxis +#' @param slow_parasite_clearance_time single value representing the mean time individual's experiencing slow parasite clearance reside in the treated state +#' @export +set_antimalarial_resistance <- function(parameters, + drug, + timesteps, + artemisinin_resistance, + partner_drug_resistance, + slow_parasite_clearance_prob, + early_treatment_failure_prob, + late_clinical_failure_prob, + late_parasitological_prob, + reinfection_prob, + slow_parasite_clearance_time) { + + if(any(c(length(artemisinin_resistance), + length(partner_drug_resistance), + length(slow_parasite_clearance_prob), + length(early_treatment_failure_prob), + length(late_clinical_failure_prob), + length(late_parasitological_prob), + length(reinfection_prob)) != length(timesteps))) { + stop("Length of one or more resistance parameter vectors does not match time steps specified for update") + } + + if(any(artemisinin_resistance < 0 | artemisinin_resistance > 1 | + partner_drug_resistance < 0 | partner_drug_resistance > 1)) { + stop("Artemisinin and partner-drug resistance proportions must fall between 0 and 1") + } + + if(any(slow_parasite_clearance_prob < 0 | slow_parasite_clearance_prob > 1 | + early_treatment_failure_prob < 0 | early_treatment_failure_prob > 1 | + late_clinical_failure_prob < 0 | late_clinical_failure_prob > 1 | + late_parasitological_prob < 0 | late_parasitological_prob > 1 | + reinfection_prob < 0 | reinfection_prob > 1)) { + stop("Resistance outcome probabilities must fall between 0 and 1") + } + + if(length(slow_parasite_clearance_time) != 1) { + stop("Error: length of slow_parasite_clearance_time not equal to 1") + } + + if(slow_parasite_clearance_time <= 0) { + stop("Error: slow_parasite_clearance_time is non-positive") + } + + parameters$antimalarial_resistance <- TRUE + + n_drugs <- length(parameters$drug_efficacy) + + if (drug < 1 | drug > n_drugs) { + stop('Drug index is invalid, please set drugs using set_drugs') + } + + drug_index <- which(parameters$antimalarial_resistance_drug == drug) + + if (length(drug_index) == 0) { + drug_index <- length(parameters$antimalarial_resistance_drug) + 1 + } + + parameters$antimalarial_resistance_drug[[drug_index]] <- drug + parameters$antimalarial_resistance_timesteps[[drug_index]] <- timesteps + parameters$prop_artemisinin_resistant[[drug_index]] <- artemisinin_resistance + parameters$prop_partner_drug_resistant[[drug_index]] <- partner_drug_resistance + parameters$slow_parasite_clearance_prob[[drug_index]] <- slow_parasite_clearance_prob + parameters$early_treatment_failure_prob[[drug_index]] <- early_treatment_failure_prob + parameters$late_clinical_failure_prob[[drug_index]] <- late_clinical_failure_prob + parameters$late_parasitological_failure_prob[[drug_index]] <- late_parasitological_prob + parameters$reinfection_during_prophylaxis[[drug_index]] <- reinfection_prob + parameters$dt_slow_parasite_clearance[[drug_index]] <- slow_parasite_clearance_time + + return(parameters) + +} + +#' @title Retrieve resistance parameters +#' @description +#' Retrieve the resistance parameters associated with the drug each individual receiving clinical +#' treatment has been administered in the current time step. +#' +#' @param parameters the model parameters +#' @param drugs vector of integers representing the drugs administered to each individual receiving treatment +#' @param timestep the current time step +get_antimalarial_resistance_parameters <- function(parameters, drugs, timestep) { + + if(!parameters$antimalarial_resistance) { + stop("Error: Antimalarial resistance has not been parameterised; antimalarial_resistance = FALSE") + } + + blank_vector <- numeric(length = length(drugs)) + artemisinin_resistance_proportion <- blank_vector + partner_drug_resistance_proportion <- blank_vector + slow_parasite_clearance_probability <- blank_vector + early_treatment_failure_probability <- blank_vector + late_clinical_failure_probability <- blank_vector + late_parasitological_failure_probability <- blank_vector + reinfection_during_prophylaxis_probability <- blank_vector + dt_slow_parasite_clearance <- rep(parameters$dt, length = length(drugs)) + + for(i in seq_along(parameters$antimalarial_resistance_drug)) { + drug <- parameters$antimalarial_resistance_drug[[i]] + treated_with_drug <- which(drugs == drug) + resistance_timestep <- match_timestep(ts = parameters$antimalarial_resistance_timesteps[[i]], t = timestep) + artemisinin_resistance_proportion[treated_with_drug] <- parameters$prop_artemisinin_resistant[[i]][resistance_timestep] + partner_drug_resistance_proportion[treated_with_drug] <- parameters$prop_partner_drug_resistant[[i]][resistance_timestep] + slow_parasite_clearance_probability[treated_with_drug] <- parameters$slow_parasite_clearance_prob[[i]][resistance_timestep] + early_treatment_failure_probability[treated_with_drug] <- parameters$early_treatment_failure_prob[[i]][resistance_timestep] + late_clinical_failure_probability[treated_with_drug] <- parameters$late_clinical_failure_prob[[i]][resistance_timestep] + late_parasitological_failure_probability[treated_with_drug] <- parameters$late_parasitological_failure_prob[[i]][resistance_timestep] + reinfection_during_prophylaxis_probability[treated_with_drug] <- parameters$reinfection_during_prophylaxis[[i]][resistance_timestep] + dt_slow_parasite_clearance[treated_with_drug] <- parameters$dt_slow_parasite_clearance[[i]] + } + + resistance_parameters <- list() + resistance_parameters$artemisinin_resistance_proportion <- artemisinin_resistance_proportion + resistance_parameters$partner_drug_resistance_proportion <- partner_drug_resistance_proportion + resistance_parameters$slow_parasite_clearance_probability <- slow_parasite_clearance_probability + resistance_parameters$early_treatment_failure_probability <- early_treatment_failure_probability + resistance_parameters$late_clinical_failure_probability <- late_clinical_failure_probability + resistance_parameters$late_parasitological_failure_probability <- late_parasitological_failure_probability + resistance_parameters$reinfection_during_prophylaxis_probability <- reinfection_during_prophylaxis_probability + resistance_parameters$dt_slow_parasite_clearance <- dt_slow_parasite_clearance + + return(resistance_parameters) + +} \ No newline at end of file diff --git a/R/human_infection.R b/R/human_infection.R index 1cb521c4..6b2ac3cf 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -263,23 +263,27 @@ update_severe_disease <- function( #' @param renderer simulation renderer #' @noRd calculate_treated <- function( - variables, - clinical_infections, - parameters, - timestep, - renderer - ) { + variables, + clinical_infections, + parameters, + timestep, + renderer +) { + + if(clinical_infections$size() == 0) { + return(individual::Bitset$new(parameters$human_population)) + } + treatment_coverages <- get_treatment_coverages(parameters, timestep) ft <- sum(treatment_coverages) - + if (ft == 0) { return(individual::Bitset$new(parameters$human_population)) } - + renderer$render('ft', ft, timestep) seek_treatment <- sample_bitset(clinical_infections, ft) n_treat <- seek_treatment$size() - renderer$render('n_treated', n_treat, timestep) drugs <- as.numeric(parameters$clinical_treatment_drugs[ @@ -290,29 +294,51 @@ calculate_treated <- function( replace = TRUE ) ]) - - successful <- bernoulli_multi_p(parameters$drug_efficacy[drugs]) - treated_index <- bitset_at(seek_treatment, successful) - - # Update those who have been treated - if (treated_index$size() > 0) { - variables$state$queue_update('Tr', treated_index) + + if(parameters$antimalarial_resistance) { + resistance_parameters <- get_antimalarial_resistance_parameters( + parameters = parameters, + drugs = drugs, + timestep = timestep + ) + unsuccessful_treatment_probability <- resistance_parameters$artemisinin_resistance_proportion * resistance_parameters$early_treatment_failure_probability + susceptible_to_treatment_index <- bernoulli_multi_p(p = 1 - unsuccessful_treatment_probability) + susceptible_to_treatment <- bitset_at(seek_treatment, susceptible_to_treatment_index) + drugs <- drugs[susceptible_to_treatment_index] + } else { + susceptible_to_treatment <- seek_treatment + + } + + n_early_treatment_failure <- n_treat - susceptible_to_treatment$size() + successfully_treated_index <- bernoulli_multi_p(parameters$drug_efficacy[drugs]) + successfully_treated <- bitset_at(susceptible_to_treatment, successfully_treated_index) + successfully_treated_drugs <- drugs[successfully_treated_index] + n_treat_eff_fail <- susceptible_to_treatment$size() - length(successfully_treated_index) + renderer$render('n_early_treatment_failure', n_early_treatment_failure, timestep) + renderer$render('n_treat_eff_fail', n_treat_eff_fail, timestep) + renderer$render('n_treat_success', successfully_treated$size(), timestep) + + # Update variables of those who have been successfully treated: + if (successfully_treated$size() > 0) { + variables$state$queue_update("Tr", successfully_treated) variables$infectivity$queue_update( - parameters$cd * parameters$drug_rel_c[drugs[successful]], - treated_index + parameters$cd * parameters$drug_rel_c[successfully_treated_drugs], + successfully_treated ) variables$drug$queue_update( - drugs[successful], - treated_index + successfully_treated_drugs, + successfully_treated ) variables$drug_time$queue_update( timestep, - treated_index + successfully_treated ) } - treated_index + successfully_treated } + #' @title Schedule infections #' @description #' Schedule infections in humans after the incubation period diff --git a/R/model.R b/R/model.R index fb6ffef9..4c4f7c2b 100644 --- a/R/model.R +++ b/R/model.R @@ -74,6 +74,9 @@ #' susceptible #' * net_usage: the number people protected by a bed net #' * mosquito_deaths: number of adult female mosquitoes who die this timestep +#' * n_early_treatment_failure: number of clinically treated individuals who experienced early treatment failure in this timestep +#' * n_treat_eff_fail: number of clinically treated individuals who's treatment failed due to drug efficacy +#' * n_treat_success: number of successfully treated individuals in this timestep #' #' @param timesteps the number of timesteps to run the simulation for (in days) #' @param parameters a named list of parameters to use diff --git a/R/parameters.R b/R/parameters.R index 8a96989a..9abbeb60 100644 --- a/R/parameters.R +++ b/R/parameters.R @@ -192,6 +192,22 @@ #' * tbv_tra_mu - transmission reduction parameter; default = 12.63 #' * tbv_gamma1 - transmission reduction parameter; default = 2.5 #' * tbv_gamma2 - transmission reduction parameter; default = 0.06 +#' +#' Antimalarial resistance parameters: +#' please set antimalarial resistance parameters with the convenience functions in +#' `antimalarial_resistance.R:set_antimalarial_resistance` +#' +#' * antimalarial_resistance - boolean for if antimalarial resistance is enabled; default = FALSE +#' * antimalarial_resistance_drug - vector of drugs for which resistance can be parameterised; default = NULL +#' * antimalarial_resistance_timesteps - vector of time steps on which resistance updates occur; default = NULL +#' * prop_artemisinin_resistant - vector of proportions of infections resistant to the artemisinin component of a given drug; default = NULL +#' * prop_partner_drug_resistant - vector of proportions of infections resistant to the parter drug component of a given drug; default = NULL +#' * slow_parasite_clearance_prob - vector of probabilities of slow parasite clearance for a given drug; default = NULL +#' * early_treatment_failure_prob - vector of probabilities of early treatment failure for a given drug; default = NULL +#' * late_clinical_failure_prob - vector of probabilities of late clinical failure for a given drug; default = NULL +#' * late_parasitological_failure_prob - vector of probabilities of late parasitological failure for a given drug; default = NULL +#' * reinfection_during_prophylaxis - vector of probabilities of reinfection during prophylaxis for a given drug; default = NULL +#' * dt_slow_parasite_clearance - the delay for humans experiencing slow parasite clearance to move from state Tr to S; default = NULL #' #' rendering: #' All values are in timesteps and all ranges are inclusive @@ -377,6 +393,17 @@ get_parameters <- function(overrides = list()) { tbv_timesteps = NULL, tbv_coverages = NULL, tbv_ages = NULL, + # antimalarial resistance + antimalarial_resistance = FALSE, + antimalarial_resistance_drug = NULL, + antimalarial_resistance_timesteps = NULL, + prop_artemisinin_resistant = NULL, + prop_partner_drug_resistant = NULL, + slow_parasite_clearance_prob = NULL, + early_treatment_failure_prob = NULL, + late_clinical_failure_prob = NULL, + late_parasitological_failure_prob = NULL, + reinfection_during_prophylaxis = NULL, # flexible carrying capacity carrying_capacity = FALSE, carrying_capacity_timesteps = NULL, diff --git a/man/CorrelationParameters.Rd b/man/CorrelationParameters.Rd index b1d22578..c2e6ada7 100644 --- a/man/CorrelationParameters.Rd +++ b/man/CorrelationParameters.Rd @@ -14,17 +14,17 @@ Describes an event in the simulation \section{Methods}{ \subsection{Public methods}{ \itemize{ -\item \href{#method-new}{\code{CorrelationParameters$new()}} -\item \href{#method-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}} -\item \href{#method-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}} -\item \href{#method-sigma}{\code{CorrelationParameters$sigma()}} -\item \href{#method-mvnorm}{\code{CorrelationParameters$mvnorm()}} -\item \href{#method-clone}{\code{CorrelationParameters$clone()}} +\item \href{#method-CorrelationParameters-new}{\code{CorrelationParameters$new()}} +\item \href{#method-CorrelationParameters-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}} +\item \href{#method-CorrelationParameters-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}} +\item \href{#method-CorrelationParameters-sigma}{\code{CorrelationParameters$sigma()}} +\item \href{#method-CorrelationParameters-mvnorm}{\code{CorrelationParameters$mvnorm()}} +\item \href{#method-CorrelationParameters-clone}{\code{CorrelationParameters$clone()}} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-new}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-CorrelationParameters-new}{}}} \subsection{Method \code{new()}}{ initialise correlation parameters \subsection{Usage}{ @@ -40,8 +40,8 @@ initialise correlation parameters } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-inter_round_rho}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_round_rho}{}}} \subsection{Method \code{inter_round_rho()}}{ Add rho between rounds \subsection{Usage}{ @@ -60,8 +60,8 @@ the intervention} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-inter_intervention_rho}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_intervention_rho}{}}} \subsection{Method \code{inter_intervention_rho()}}{ Add rho between interventions \subsection{Usage}{ @@ -83,8 +83,8 @@ the intervention} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-sigma}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-CorrelationParameters-sigma}{}}} \subsection{Method \code{sigma()}}{ Standard deviation of each intervention between rounds \subsection{Usage}{ @@ -93,8 +93,8 @@ Standard deviation of each intervention between rounds } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-mvnorm}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-CorrelationParameters-mvnorm}{}}} \subsection{Method \code{mvnorm()}}{ multivariate norm draws for these parameters \subsection{Usage}{ @@ -103,8 +103,8 @@ multivariate norm draws for these parameters } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-clone}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-CorrelationParameters-clone}{}}} \subsection{Method \code{clone()}}{ The objects of this class are cloneable with this method. \subsection{Usage}{ diff --git a/man/get_antimalarial_resistance_parameters.Rd b/man/get_antimalarial_resistance_parameters.Rd new file mode 100644 index 00000000..71de6961 --- /dev/null +++ b/man/get_antimalarial_resistance_parameters.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/antimalarial_resistance.R +\name{get_antimalarial_resistance_parameters} +\alias{get_antimalarial_resistance_parameters} +\title{Retrieve resistance parameters} +\usage{ +get_antimalarial_resistance_parameters(parameters, drugs, timestep) +} +\arguments{ +\item{parameters}{the model parameters} + +\item{drugs}{vector of integers representing the drugs administered to each individual receiving treatment} + +\item{timestep}{the current time step} +} +\description{ +Retrieve the resistance parameters associated with the drug each individual receiving clinical +treatment has been administered in the current time step. +} diff --git a/man/get_parameters.Rd b/man/get_parameters.Rd index 54188c01..fc6ba000 100644 --- a/man/get_parameters.Rd +++ b/man/get_parameters.Rd @@ -214,6 +214,23 @@ please set TBV parameters with the convenience functions in \item tbv_gamma2 - transmission reduction parameter; default = 0.06 } +Antimalarial resistance parameters: +please set antimalarial resistance parameters with the convenience functions in +\code{antimalarial_resistance.R:set_antimalarial_resistance} +\itemize{ +\item antimalarial_resistance - boolean for if antimalarial resistance is enabled; default = FALSE +\item antimalarial_resistance_drug - vector of drugs for which resistance can be parameterised; default = NULL +\item antimalarial_resistance_timesteps - vector of time steps on which resistance updates occur; default = NULL +\item prop_artemisinin_resistant - vector of proportions of infections resistant to the artemisinin component of a given drug; default = NULL +\item prop_partner_drug_resistant - vector of proportions of infections resistant to the parter drug component of a given drug; default = NULL +\item slow_parasite_clearance_prob - vector of probabilities of slow parasite clearance for a given drug; default = NULL +\item early_treatment_failure_prob - vector of probabilities of early treatment failure for a given drug; default = NULL +\item late_clinical_failure_prob - vector of probabilities of late clinical failure for a given drug; default = NULL +\item late_parasitological_failure_prob - vector of probabilities of late parasitological failure for a given drug; default = NULL +\item reinfection_during_prophylaxis - vector of probabilities of reinfection during prophylaxis for a given drug; default = NULL +\item dt_slow_parasite_clearance - the delay for humans experiencing slow parasite clearance to move from state Tr to S; default = NULL +} + rendering: All values are in timesteps and all ranges are inclusive \itemize{ diff --git a/man/run_simulation.Rd b/man/run_simulation.Rd index 64cac9c4..07871b13 100644 --- a/man/run_simulation.Rd +++ b/man/run_simulation.Rd @@ -90,5 +90,8 @@ subpatent susceptible \item net_usage: the number people protected by a bed net \item mosquito_deaths: number of adult female mosquitoes who die this timestep +\item n_early_treatment_failure: number of clinically treated individuals who experienced early treatment failure in this timestep +\item n_treat_eff_fail: number of clinically treated individuals who's treatment failed due to drug efficacy +\item n_treat_success: number of successfully treated individuals in this timestep } } diff --git a/man/set_antimalarial_resistance.Rd b/man/set_antimalarial_resistance.Rd new file mode 100644 index 00000000..82c65405 --- /dev/null +++ b/man/set_antimalarial_resistance.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/antimalarial_resistance.R +\name{set_antimalarial_resistance} +\alias{set_antimalarial_resistance} +\title{Parameterise antimalarial resistance} +\usage{ +set_antimalarial_resistance( + parameters, + drug, + timesteps, + artemisinin_resistance, + partner_drug_resistance, + slow_parasite_clearance_prob, + early_treatment_failure_prob, + late_clinical_failure_prob, + late_parasitological_prob, + reinfection_prob, + slow_parasite_clearance_time +) +} +\arguments{ +\item{parameters}{the model parameters} + +\item{drug}{the index of the drug which resistance is being set, as set by the set_drugs() function, in the parameter list} + +\item{timesteps}{vector of time steps for each update to resistance proportion and resistance outcome probability} + +\item{artemisinin_resistance}{vector of updates to the proportions of infections that are artemisinin resistant at time t} + +\item{partner_drug_resistance}{vector of updates to the proportions of infections that are partner-drug resistant at time t} + +\item{slow_parasite_clearance_prob}{vector of updates to the proportion of artemisinin-resistant infections that result in early treatment failure} + +\item{early_treatment_failure_prob}{vector of updates to the proportion of artemisinin-resistant infections that result in slow parasite clearance} + +\item{late_clinical_failure_prob}{vector of updates to the proportion of partner-drug-resistant infections that result in late clinical failure} + +\item{late_parasitological_prob}{vector of updates to the proportion of partner-drug-resistant infections that result in late parasitological failure} + +\item{reinfection_prob}{vector of updates to the proportion of partner-drug-resistant infections that result in reinfection during prophylaxis} + +\item{slow_parasite_clearance_time}{single value representing the mean time individual's experiencing slow parasite clearance reside in the treated state} +} +\description{ +Parameterise antimalarial resistance +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index affb233d..f5c226fd 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -260,7 +260,7 @@ BEGIN_RCPP END_RCPP } -RcppExport SEXP run_testthat_tests(void); +RcppExport SEXP run_testthat_tests(); static const R_CallMethodDef CallEntries[] = { {"_malariasimulation_create_adult_mosquito_model", (DL_FUNC) &_malariasimulation_create_adult_mosquito_model, 5}, diff --git a/tests/testthat/test-antimalarial-resistance.R b/tests/testthat/test-antimalarial-resistance.R new file mode 100644 index 00000000..400d93a3 --- /dev/null +++ b/tests/testthat/test-antimalarial-resistance.R @@ -0,0 +1,322 @@ +test_that('set_antimalarial_resistance() toggles resistance on', { + simparams <- get_parameters() + simparams <- set_drugs(parameters = simparams, drugs = list(SP_AQ_params)) + simparams <- set_clinical_treatment(parameters = simparams, + drug = 1, + timesteps = 1, + coverages = 1) + set_antimalarial_resistance(parameters = simparams, + drug = 1, + timesteps = 1, + artemisinin_resistance = 0.5, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0.5, + early_treatment_failure_prob = 0.6, + late_clinical_failure_prob = 0.2, + late_parasitological_prob = 0.3, + reinfection_prob = 0.4, + slow_parasite_clearance_time = 10) -> simparams + expect_identical(object = simparams$antimalarial_resistance, expected = TRUE) +}) + +test_that('set_antimalarial_resistance() errors if parameter inputs of different length to timesteps', { + simparams <- get_parameters() + simparams <- set_drugs(parameters = simparams, drugs = list(SP_AQ_params)) + simparams <- set_clinical_treatment(parameters = simparams, + drug = 1, + timesteps = 1, + coverages = 1) + expect_error(object = set_antimalarial_resistance(parameters = simparams, + drug = 1, + timesteps = c(1, 10), + artemisinin_resistance = 0.5, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0.5, + early_treatment_failure_prob = 0.6, + late_clinical_failure_prob = 0.2, + late_parasitological_prob = 0.3, + reinfection_prob = 0.4, + slow_parasite_clearance_time = 10)) +}) + +test_that('set_antimalarial_resistance() errors if resistance proportions outside of range 0-1', { + simparams <- get_parameters() + simparams <- set_drugs(parameters = simparams, drugs = list(SP_AQ_params)) + simparams <- set_clinical_treatment(parameters = simparams, + drug = 1, + timesteps = 1, + coverages = 1) + expect_error(object = set_antimalarial_resistance(parameters = simparams, + drug = 1, + timesteps = 1, + artemisinin_resistance = 1.01, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0.5, + early_treatment_failure_prob = 0.6, + late_clinical_failure_prob = 0.2, + late_parasitological_prob = 0.3, + reinfection_prob = 0.4, + slow_parasite_clearance_time = 10), + regexp = "Artemisinin and partner-drug resistance proportions must fall between 0 and 1") +}) + +test_that('set_antimalarial_resistance() errors if resistance phenotype probabilities outside bound of 0-1', { + simparams <- get_parameters() + simparams <- set_drugs(parameters = simparams, drugs = list(SP_AQ_params)) + simparams <- set_clinical_treatment(parameters = simparams, + drug = 1, + timesteps = 1, + coverages = 1) + expect_error(object = set_antimalarial_resistance(parameters = simparams, + drug = 1, + timesteps = 1, + artemisinin_resistance = 0.4, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = -0.5, + early_treatment_failure_prob = 0.6, + late_clinical_failure_prob = 0.2, + late_parasitological_prob = 0.3, + reinfection_prob = 0.4, + slow_parasite_clearance_time = 5)) +}) + +test_that('set_antimalarial_resistance() errors if drug index > than number of drugs assigned using set_drugs()', { + simparams <- get_parameters() + simparams <- set_drugs(parameters = simparams, drugs = list(SP_AQ_params)) + simparams <- set_clinical_treatment(parameters = simparams, + drug = 1, + timesteps = 1, + coverages = 1) + expect_error(object = set_antimalarial_resistance(parameters = simparams, + drug = 2, + timesteps = 1, + artemisinin_resistance = 0.4, + partner_drug_resistance = 0.3, + slow_parasite_clearance_prob = 0.5, + early_treatment_failure_prob = 0.6, + late_clinical_failure_prob = 0.2, + late_parasitological_prob = 0.3, + reinfection_prob = 0.4)) +}) + +test_that('set_antimalarial_resistance() assigns parameters correctly despite order in which resistance parameters are specified', { + + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(AL_params, SP_AQ_params, DHA_PQP_params)) + parameters <- set_clinical_treatment(parameters = parameters, drug = 2, timesteps = 1, coverages = 0.2) + parameters <- set_clinical_treatment(parameters = parameters, drug = 1, timesteps = 1, coverages = 0.1) + parameters <- set_clinical_treatment(parameters = parameters, drug = 3, timesteps = 1, coverages = 0.4) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 2, + timesteps = 1, + artemisinin_resistance = 0.5, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0.41, + early_treatment_failure_prob = 0.2, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 5) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 3, + timesteps = 1, + artemisinin_resistance = 0, + partner_drug_resistance = 0.43, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 0, + late_clinical_failure_prob = 0.01, + late_parasitological_prob = 0.42, + reinfection_prob = 0.89, + slow_parasite_clearance_time = 10) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 1, + timesteps = 1, + artemisinin_resistance = 0.27, + partner_drug_resistance = 0.61, + slow_parasite_clearance_prob = 0.23, + early_treatment_failure_prob = 0.9, + late_clinical_failure_prob = 0.49, + late_parasitological_prob = 0.81, + reinfection_prob = 0.009, + slow_parasite_clearance_time = 20) + + expect_identical(parameters$antimalarial_resistance, TRUE) + expect_identical(unlist(parameters$antimalarial_resistance_drug), c(2, 3, 1)) + expect_identical(unlist(parameters$antimalarial_resistance_timesteps), rep(1, 3)) + expect_identical(unlist(parameters$prop_artemisinin_resistant), c(0.5, 0, 0.27)) + expect_identical(unlist(parameters$prop_partner_drug_resistant), c(0, 0.43, 0.61)) + expect_identical(unlist(parameters$slow_parasite_clearance_prob), c(0.41, 0, 0.23)) + expect_identical(unlist(parameters$early_treatment_failure_prob), c(0.2, 0, 0.9)) + expect_identical(unlist(parameters$late_clinical_failure_prob), c(0, 0.01, 0.49)) + expect_identical(unlist(parameters$late_parasitological_failure_prob), c(0, 0.42, 0.81)) + expect_identical(unlist(parameters$reinfection_during_prophylaxis), c(0, 0.89, 0.009)) + expect_identical(unlist(parameters$dt_slow_parasite_clearance), c(5, 10, 20)) + +}) + +test_that(desc = "set_antimalarial_resistance errors if length slow_parasite_clearance_time > 1", code = { + + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(SP_AQ_params)) + parameters <- set_clinical_treatment(parameters = parameters, + drug = 1, + timesteps = c(0, 10), + coverages = c(0.1, 0.2)) + + expect_error( + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 1, + timesteps = c(0, 10), + artemisinin_resistance = c(0.4, 0.8), + partner_drug_resistance = c(0.23, 0.43), + slow_parasite_clearance_prob = c(0.2, 0.4), + early_treatment_failure_prob = c(0, 0.45), + late_clinical_failure_prob = c(0.01, 0.01), + late_parasitological_prob = c(0.05, 0.06), + reinfection_prob = c(0.86, 0.86), + slow_parasite_clearance_time = c(10 ,11)), + "Error: length of slow_parasite_clearance_time not equal to 1") +}) + +test_that(desc = "set_antimalarial_resistance errors if slow_parasite_clearance_time not positive", code = { + + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(SP_AQ_params)) + parameters <- set_clinical_treatment(parameters = parameters, + drug = 1, + timesteps = c(0, 10), + coverages = c(0.1, 0.2)) + + expect_error( + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 1, + timesteps = c(0, 10), + artemisinin_resistance = c(0.4, 0.8), + partner_drug_resistance = c(0.23, 0.43), + slow_parasite_clearance_prob = c(0.2, 0.4), + early_treatment_failure_prob = c(0, 0.45), + late_clinical_failure_prob = c(0.01, 0.01), + late_parasitological_prob = c(0.05, 0.06), + reinfection_prob = c(0.86, 0.86), + slow_parasite_clearance_time = c(0)), + "Error: slow_parasite_clearance_time is non-positive") +}) + +test_that('get_antimalarial_resistance_parameters() correctly retrieves parameters when multiple drugs assigned', { + + get_parameters(overrides = list(human_population = 10000)) %>% + set_drugs(drugs = list(AL_params, SP_AQ_params, DHA_PQP_params)) %>% + set_clinical_treatment(drug = 1, timesteps = 1, coverages = 0.4) %>% + set_clinical_treatment(drug = 2, timesteps = 1, coverages = 0.3) %>% + set_clinical_treatment(drug = 3, timesteps = 1, coverages = 0.2) %>% + set_equilibrium(init_EIR = 20) %>% + set_antimalarial_resistance(drug = 2, + timesteps = c(0, 20), + artemisinin_resistance = c(0.02, 0.2), + partner_drug_resistance = c(0.02, 0.2), + slow_parasite_clearance_prob = c(0.02, 0.2), + early_treatment_failure_prob = c(0.02, 0.2), + late_clinical_failure_prob = c(0.02, 0.2), + late_parasitological_prob = c(0.02, 0.2), + reinfection_prob = c(0.02, 0.2), + slow_parasite_clearance_time = 20) %>% + set_antimalarial_resistance(drug = 1, + timesteps = c(0, 10), + artemisinin_resistance = c(0.01, 0.1), + partner_drug_resistance = c(0.01, 0.1), + slow_parasite_clearance_prob = c(0.01, 0.1), + early_treatment_failure_prob = c(0.01, 0.1), + late_clinical_failure_prob = c(0.01, 0.1), + late_parasitological_prob = c(0.01, 0.1), + reinfection_prob = c(0.01, 0.1), + slow_parasite_clearance_time = 10) %>% + set_antimalarial_resistance(drug = 3, + timesteps = c(0, 30), + artemisinin_resistance = c(0.03, 0.3), + partner_drug_resistance = c(0.03, 0.3), + slow_parasite_clearance_prob = c(0.03, 0.3), + early_treatment_failure_prob = c(0.03, 0.3), + late_clinical_failure_prob = c(0.03, 0.3), + late_parasitological_prob = c(0.03, 0.3), + reinfection_prob = c(0.03, 0.3), + slow_parasite_clearance_time = 30) -> parameters + + drugs <- c(1, 3, 2, 1, 2, 3, 3, 3, 2, 1, 3, 1, 2, 3, 2) + timestep <- 25 + + resistance_parameters <- get_antimalarial_resistance_parameters(parameters = parameters, + drugs = drugs, + timestep = timestep) + + expected_resistance_parameters <- list() + expected_resistance_parameters$artemisinin_resistance_proportion <- c(0.1, 0.03, 0.2, 0.1, 0.2, 0.03, 0.03, 0.03, 0.2, 0.1, 0.03, 0.1, 0.2, 0.03, 0.2) + expected_resistance_parameters$partner_drug_resistance_proportion <- c(0.1, 0.03, 0.2, 0.1, 0.2, 0.03, 0.03, 0.03, 0.2, 0.1, 0.03, 0.1, 0.2, 0.03, 0.2) + expected_resistance_parameters$slow_parasite_clearance_probability <- c(0.1, 0.03, 0.2, 0.1, 0.2, 0.03, 0.03, 0.03, 0.2, 0.1, 0.03, 0.1, 0.2, 0.03, 0.2) + expected_resistance_parameters$early_treatment_failure_probability <- c(0.1, 0.03, 0.2, 0.1, 0.2, 0.03, 0.03, 0.03, 0.2, 0.1, 0.03, 0.1, 0.2, 0.03, 0.2) + expected_resistance_parameters$late_clinical_failure_probability <- c(0.1, 0.03, 0.2, 0.1, 0.2, 0.03, 0.03, 0.03, 0.2, 0.1, 0.03, 0.1, 0.2, 0.03, 0.2) + expected_resistance_parameters$late_parasitological_failure_probability <- c(0.1, 0.03, 0.2, 0.1, 0.2, 0.03, 0.03, 0.03, 0.2, 0.1, 0.03, 0.1, 0.2, 0.03, 0.2) + expected_resistance_parameters$reinfection_during_prophylaxis_probability <- c(0.1, 0.03, 0.2, 0.1, 0.2, 0.03, 0.03, 0.03, 0.2, 0.1, 0.03, 0.1, 0.2, 0.03, 0.2) + expected_resistance_parameters$dt_slow_parasite_clearance <- c(10, 30, 20, 10, 20, 30, 30, 30, 20, 10, 30, 10, 20, 30, 20) + + expect_identical(resistance_parameters, expected = expected_resistance_parameters) + + }) + +test_that('get_antimalarial_resistance_parameters() correctly retrieves parameters when not all drugs assigned resistance', { + + get_parameters(overrides = list(human_population = 10000)) %>% + set_drugs(drugs = list(AL_params, SP_AQ_params, DHA_PQP_params)) %>% + set_clinical_treatment(drug = 1, timesteps = 1, coverages = 0.4) %>% + set_clinical_treatment(drug = 2, timesteps = 1, coverages = 0.3) %>% + set_clinical_treatment(drug = 3, timesteps = 1, coverages = 0.2) %>% + set_equilibrium(init_EIR = 20) %>% + set_antimalarial_resistance(drug = 2, + timesteps = c(0, 20), + artemisinin_resistance = c(0.02, 0.2), + partner_drug_resistance = c(0.02, 0.2), + slow_parasite_clearance_prob = c(0.02, 0.2), + early_treatment_failure_prob = c(0.02, 0.2), + late_clinical_failure_prob = c(0.02, 0.2), + late_parasitological_prob = c(0.02, 0.2), + reinfection_prob = c(0.02, 0.2), + slow_parasite_clearance_time = 20) -> parameters + + drugs <- c(1, 3, 2, 1, 2, 3, 3, 3, 2, 1, 3, 1, 2, 3, 2) + timestep <- 25 + + resistance_parameters <- get_antimalarial_resistance_parameters(parameters = parameters, + drugs = drugs, + timestep = timestep) + + expected_resistance_parameters <- list() + expected_resistance_parameters$artemisinin_resistance_proportion <- c(0, 0, 0.2, 0, 0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2, 0, 0.2) + expected_resistance_parameters$partner_drug_resistance_proportion <- c(0, 0, 0.2, 0, 0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2, 0, 0.2) + expected_resistance_parameters$slow_parasite_clearance_probability <- c(0, 0, 0.2, 0, 0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2, 0, 0.2) + expected_resistance_parameters$early_treatment_failure_probability <- c(0, 0, 0.2, 0, 0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2, 0, 0.2) + expected_resistance_parameters$late_clinical_failure_probability <- c(0, 0, 0.2, 0, 0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2, 0, 0.2) + expected_resistance_parameters$late_parasitological_failure_probability <- c(0, 0, 0.2, 0, 0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2, 0, 0.2) + expected_resistance_parameters$reinfection_during_prophylaxis_probability <- c(0, 0, 0.2, 0, 0.2, 0, 0, 0, 0.2, 0, 0, 0, 0.2, 0, 0.2) + expected_resistance_parameters$dt_slow_parasite_clearance <- c(5, 5, 20, 5, 20, 5, 5, 5, 20, 5, 5, 5, 20, 5, 20) + + expect_identical(resistance_parameters, expected = expected_resistance_parameters) + + }) + +test_that('get_antimalarial_resistance_parameters() returns an error when antimalarial resistance has not been parameterised', { + + get_parameters(overrides = list(human_population = 10000)) %>% + set_drugs(drugs = list(AL_params, SP_AQ_params, DHA_PQP_params)) %>% + set_clinical_treatment(drug = 1, timesteps = 1, coverages = 0.4) %>% + set_clinical_treatment(drug = 2, timesteps = 1, coverages = 0.3) %>% + set_clinical_treatment(drug = 3, timesteps = 1, coverages = 0.2) %>% + set_equilibrium(init_EIR = 20) -> parameters + + drugs <- c(1, 3, 2, 1, 2, 3, 3, 3, 2, 1, 3, 1, 2, 3, 2) + timestep <- 25 + + + expect_error(get_antimalarial_resistance_parameters(parameters = parameters, + drugs = drugs, + timestep = timestep), + "Error: Antimalarial resistance has not been parameterised; antimalarial_resistance = FALSE") +}) \ No newline at end of file diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index 0400a361..70a0ec6f 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -206,7 +206,6 @@ test_that('calculate_infections works various combinations of drug and vaccinati }) - test_that('calculate_clinical_infections correctly samples clinically infected', { timestep <- 5 parameters <- get_parameters() @@ -265,10 +264,10 @@ test_that('calculate_treated correctly samples treated and updates the drug stat drug = list(queue_update = mockery::mock()), drug_time = list(queue_update = mockery::mock()) ) - + recovery_mock <- mockery::mock() mockery::stub(calculate_treated, 'recovery$schedule', recovery_mock) - + seek_treatment <- individual::Bitset$new(4)$insert(c(1, 2, 4)) mockery::stub( calculate_treated, @@ -280,7 +279,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat bernoulli_mock <- mockery::mock(c(1, 3)) mockery::stub(calculate_treated, 'bernoulli_multi_p', bernoulli_mock) mockery::stub(calculate_treated, 'log_uniform', mockery::mock(c(3, 4))) - + clinical_infections <- individual::Bitset$new(4) clinical_infections$insert(c(1, 2, 3, 4)) calculate_treated( @@ -290,7 +289,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat timestep, mock_render(timestep) ) - + mockery::expect_args( sample_mock, 1, @@ -299,13 +298,13 @@ test_that('calculate_treated correctly samples treated and updates the drug stat c(.25, .25), TRUE ) - + mockery::expect_args( bernoulli_mock, 1, parameters$drug_efficacy[c(2, 1, 1, 1)] ) - + expect_bitset_update(variables$state$queue_update, 'Tr', c(1, 4)) expect_bitset_update( variables$infectivity$queue_update, @@ -316,6 +315,236 @@ test_that('calculate_treated correctly samples treated and updates the drug stat expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 4)) }) +test_that('calculate_treated correctly samples treated and updates the drug state when resistance set', { + + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(AL_params, SP_AQ_params)) + parameters <- set_clinical_treatment(parameters = parameters, drug = 1, timesteps = 1, coverages = 0.25) + parameters <- set_clinical_treatment(parameters = parameters, drug = 2, timesteps = 1, coverages = 0.25) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 1, + timesteps = 1, + artemisinin_resistance = 0.5, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 0.2, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 10) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 2, + timesteps = 1, + artemisinin_resistance = 0, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 0.9, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 15) + + clinical_infections <- individual::Bitset$new(20)$insert(1:20) + timestep <- 5 + events <- create_events(parameters) + variables <- list( + state = list(queue_update = mockery::mock()), + infectivity = list(queue_update = mockery::mock()), + drug = list(queue_update = mockery::mock()), + drug_time = list(queue_update = mockery::mock()) + ) + + seek_treatment <- individual::Bitset$new(20)$insert(c(1:10)) + seek_treatment_mock <- mockery::mock(seek_treatment) + mockery::stub(where = calculate_treated, what = 'sample_bitset', how = seek_treatment_mock) + + mock_drugs <- mockery::mock(c(2, 1, 1, 1, 2, 2, 2, 1, 2, 1)) + mockery::stub(calculate_treated, 'sample.int', mock_drugs) + + bernoulli_mock <- mockery::mock(c(1, 2, 3, 4, 5, 6, 7, 8, 9), c(1, 2, 3, 4, 5, 6, 7)) + mockery::stub(calculate_treated, 'bernoulli_multi_p', bernoulli_mock) + + calculate_treated( + variables, + clinical_infections, + parameters, + timestep, + mock_render(timestep) + ) + + mockery::expect_args( + mock_drugs, + 1, + 2, + 10, + c(.25, .25), + TRUE + ) + + mockery::expect_args( + seek_treatment_mock, + 1, + clinical_infections, + 0.5 + ) + + mockery::expect_args( + bernoulli_mock, + 1, + c(1.0, 0.9, 0.9, 0.9, 1.0, 1.0, 1.0, 0.9, 1, 0.9) + ) + + mockery::expect_args( + bernoulli_mock, + 2, + parameters$drug_efficacy[c(2, 1, 1, 1, 2, 2, 2, 1, 2)] + ) + + expect_bitset_update(variables$state$queue_update, 'Tr', c(1, 2, 3, 4, 5, 6, 7)) + expect_bitset_update( + variables$infectivity$queue_update, + parameters$cd * parameters$drug_rel_c[c(2, 1, 1, 1, 2, 2, 2)], + c(1, 2, 3, 4, 5, 6, 7) + ) + expect_bitset_update(variables$drug$queue_update, c(2, 1, 1, 1, 2, 2, 2), c(1, 2, 3, 4, 5, 6, 7)) + expect_bitset_update(variables$drug_time$queue_update, 5, c(1, 2, 3, 4, 5, 6, 7)) + +}) + +test_that('calculate_treated correctly samples treated and updates the drug state when resistance not set for all drugs', { + + # Establish the parameters + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(AL_params, SP_AQ_params)) + parameters <- set_clinical_treatment(parameters = parameters, drug = 1, timesteps = 1, coverages = 0.25) + parameters <- set_clinical_treatment(parameters = parameters, drug = 2, timesteps = 1, coverages = 0.25) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 2, + timesteps = 1, + artemisinin_resistance = 0.8, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 0.9, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 20) + + # Establish Bitset of clinically infected individuals + clinical_infections <- individual::Bitset$new(20)$insert(1:20) + + # Set the timestep to 5: + timestep <- 5 + + # Establish the events: + events <- create_events(parameters) + + # Establish list of variables used in calculate_treated() using mocks: + variables <- list( + state = list(queue_update = mockery::mock()), + infectivity = list(queue_update = mockery::mock()), + drug = list(queue_update = mockery::mock()), + drug_time = list(queue_update = mockery::mock()) + ) + + # Create a Bitset of individuals seeking treatment individuals: + seek_treatment <- individual::Bitset$new(20)$insert(c(1:10)) + + # Create a mock of seek_treatment: + seek_treatment_mock <- mockery::mock(seek_treatment) + + # Specify that, when calculate_treated() calls sample_bitset(), return the seek_treatment_mock: + mockery::stub(where = calculate_treated, what = 'sample_bitset', how = seek_treatment_mock) + + # Create a mock_drugs object (5 of each drug): + mock_drugs <- mockery::mock(c(2, 1, 1, 1, 2, 2, 2, 1, 2, 1)) + + # Specify that when calculate_treated() calls sample.int(), it returns mock_drugs: + mockery::stub(calculate_treated, 'sample.int', mock_drugs) + + # Create a bernoulli_mock of i) individuals susceptible, and ii) individuals successfully treated: + bernoulli_mock <- mockery::mock(c(1, 2, 3, 4, 5, 6, 7, 8, 9), c(1, 2, 3, 4, 5, 6, 7)) + + # Specify that when calculate_treated() calls bernoulli_multi_p() it returns the bernoulli_mock: + mockery::stub(calculate_treated, 'bernoulli_multi_p', bernoulli_mock) + + # Run the calculate_treated() function now the mocks and stubs are established: + calculate_treated( + variables, + clinical_infections, + parameters, + timestep, + mock_render(timestep) + ) + + # Check that mock_drugs was called only once, and that the arguments used in the function call + # mock_drugs() was used in (sample.int()) match those expected: + mockery::expect_args( + mock_drugs, + 1, + 2, + 10, + c(.25, .25), + TRUE + ) + + # Check that seek_treatment_mock was called only once, and that the arguments used in the function + # call mock_drugs() was used in (sample_bitset()) match those expected: + mockery::expect_args( + seek_treatment_mock, + 1, + clinical_infections, + 0.5 + ) + + # Check that the first time bernoulli_mock was called the arguments used in the function + # call bernoulli_mock was involved in (bernoulli_multi_p()) match those expected: + mockery::expect_args( + bernoulli_mock, + 1, + c(0.28, 1, 1, 1, 0.28, 0.28, 0.28, 1, 0.28, 1) # (1 - (art_prop * etf_prob)) + ) + + # Check that the secnd time bernoulli_mock was called (bernoulli_multi_p()) the arguments used in + # the function it was called in are as expected: + mockery::expect_args( + bernoulli_mock, + 2, + parameters$drug_efficacy[c(2, 1, 1, 1, 2, 2, 2, 1, 2)] + ) + + # Check that update queued that updates the state of successfully treated individuals to "Tr" + expect_bitset_update( + variables$state$queue_update, + 'Tr', + c(1, 2, 3, 4, 5, 6, 7) + ) + + # Check that update queued that updates the infectivity of successfully treated individuals to "Tr" + # to their new infectivity (drug concentration x infectivity of "D" compartment) + expect_bitset_update( + variables$infectivity$queue_update, + parameters$cd * parameters$drug_rel_c[c(2, 1, 1, 1, 2, 2, 2)], + c(1, 2, 3, 4, 5, 6, 7) + ) + + # Check that update queued that updates the drug of successfully treated individuals to the drug + # they took: + expect_bitset_update( + variables$drug$queue_update, + c(2, 1, 1, 1, 2, 2, 2), + c(1, 2, 3, 4, 5, 6, 7) + ) + + # Check that update queued that updates the drug time of successfully treated individuals to the + # simulated/mocked time step (5) + expect_bitset_update( + variables$drug_time$queue_update, + 5, + c(1, 2, 3, 4, 5, 6, 7) + ) +}) + test_that('schedule_infections correctly schedules new infections', { parameters <- get_parameters(list(human_population = 20)) variables <- create_variables(parameters) @@ -531,3 +760,193 @@ test_that('update_severe_disease renders with no infections', { timestep ) }) + +test_that('calculate_treated returns empty Bitset when there is no clinical treatment coverage', { + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(AL_params)) + parameters <- set_clinical_treatment(parameters = parameters, drug = 1, timesteps = 1, coverages = 0) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 1, + timesteps = 1, + artemisinin_resistance = 0.5, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 0.2, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 10) + clinical_infections <- individual::Bitset$new(20)$insert(1:20) + timestep <- 5 + events <- create_events(parameters) + variables <- list( + state = list(queue_update = mockery::mock()), + infectivity = list(queue_update = mockery::mock()), + drug = list(queue_update = mockery::mock()), + drug_time = list(queue_update = mockery::mock()) + ) + renderer <- individual::Render$new(timesteps = 10) + + treated <- calculate_treated(variables = variables, + clinical_infections = clinical_infections, + parameters = parameters, + timestep = timestep, + renderer = renderer) + + expect_identical(object = treated$size(), expected = 0, info = "Error: calculate_treated() returning non-zero number of treated individuals + in the absence of clinical treatment") +}) + +test_that('calculate_treated returns empty Bitset when the clinically_infected input is an empty Bitset', { + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(AL_params)) + parameters <- set_clinical_treatment(parameters = parameters, drug = 1, timesteps = 1, coverages = 1) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 1, + timesteps = 1, + artemisinin_resistance = 0.5, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 0.2, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 10) + clinical_infections <- individual::Bitset$new(20) + timestep <- 5 + events <- create_events(parameters) + variables <- list( + state = list(queue_update = mockery::mock()), + infectivity = list(queue_update = mockery::mock()), + drug = list(queue_update = mockery::mock()), + drug_time = list(queue_update = mockery::mock()) + ) + renderer <- individual::Render$new(timesteps = 10) + + treated <- calculate_treated(variables = variables, + clinical_infections = clinical_infections, + parameters = parameters, + timestep = timestep, + renderer = renderer) + + expect_identical(object = treated$size(), expected = 0, info = "Error: calculate_treated() returning non-zero number of treated individuals + in the absence of clinically infected individuals") +}) + +test_that('calculate_treated() returns an empty Bitset when the parameter list contains no clinical + treatment or resistance parameters', { + parameters <- get_parameters() + clinical_infections <- individual::Bitset$new(20)$insert(1:20) + timestep <- 5 + events <- create_events(parameters) + variables <- list( + state = list(queue_update = mockery::mock()), + infectivity = list(queue_update = mockery::mock()), + drug = list(queue_update = mockery::mock()), + drug_time = list(queue_update = mockery::mock()) + ) + renderer <- individual::Render$new(timesteps = 10) + + treated <- calculate_treated(variables = variables, + clinical_infections = clinical_infections, + parameters = parameters, + timestep = timestep, + renderer = renderer) + + expect_identical(object = treated$size(), expected = 0, info = "Error: calculate_treated() returning non-zero number of treated individuals + in the absence of clinical treatment or resistance parameters") +}) + +test_that('Number of treatment failures matches number of individuals treated when artemisinin resistance proportion and + early treatment failure probability both set to 1', { + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(AL_params, SP_AQ_params)) + parameters <- set_clinical_treatment(parameters = parameters, + drug = 1, + timesteps = 1, + coverages = round(runif(1, 0, 1/2), + digits = 2)) + parameters <- set_clinical_treatment(parameters = parameters, + drug = 2, + timesteps = 1, + coverages = round(runif(1, 0, 1/2), + digits = 2)) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 1, + timesteps = 1, + artemisinin_resistance = 1, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 1, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 10) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 2, + timesteps = 1, + artemisinin_resistance = 1, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 1, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 20) + + clinical_infections <- individual::Bitset$new(100) + clinical_infections$insert(sample.int(n = 100, size = round(runif(n = 1, min = 10, max = 100)), replace = FALSE)) + timestep <- 5 + events <- create_events(parameters) + variables <- create_variables(parameters = parameters) + renderer <- individual::Render$new(timesteps = 10) + + treated <- calculate_treated(variables = variables, + clinical_infections = clinical_infections, + parameters = parameters, + timestep = timestep, + renderer = renderer) + + expect_identical(renderer$to_dataframe()[timestep,'n_early_treatment_failure'], renderer$to_dataframe()[timestep,'n_treated'], info = "Error: Number of + early treatment failures does not match number of treated individuals when artemisinin resistance proportion and + and early treatment failure probability both equal 1") +}) + +test_that('calculate_treated() successfully adds additional resistance columns to the renderer', { + parameters <- get_parameters() + parameters <- set_drugs(parameters = parameters, drugs = list(AL_params)) + parameters <- set_clinical_treatment(parameters = parameters, drug = 1, timesteps = 1, coverages = 1) + parameters <- set_antimalarial_resistance(parameters = parameters, + drug = 1, + timesteps = 1, + artemisinin_resistance = 0.5, + partner_drug_resistance = 0, + slow_parasite_clearance_prob = 0, + early_treatment_failure_prob = 0.5, + late_clinical_failure_prob = 0, + late_parasitological_prob = 0, + reinfection_prob = 0, + slow_parasite_clearance_time = 10) + + clinical_infections <- individual::Bitset$new(20)$insert(1:20) + timestep <- 5 + events <- create_events(parameters) + variables <- list( + state = list(queue_update = mockery::mock()), + infectivity = list(queue_update = mockery::mock()), + drug = list(queue_update = mockery::mock()), + drug_time = list(queue_update = mockery::mock()) + ) + renderer <- individual::Render$new(timesteps = 10) + + treated <- calculate_treated(variables = variables, + clinical_infections = clinical_infections, + parameters = parameters, + timestep = timestep, + renderer = renderer) + + calculate_treated_column_names <- c("ft", "n_treated", "n_early_treatment_failure", "n_treat_eff_fail", "n_treat_success") + expect_identical(sum(calculate_treated_column_names %in% colnames(renderer$to_dataframe())), length(calculate_treated_column_names), + "calculate_treated() not renderering all resistance columns when resistance is present, clinical treatment coverage + is non-zero, and the Bitset of clinically_infected individuals input is of non-zero length.") +}) diff --git a/vignettes/Antimalarial_Resistance.Rmd b/vignettes/Antimalarial_Resistance.Rmd new file mode 100644 index 00000000..7bf197b7 --- /dev/null +++ b/vignettes/Antimalarial_Resistance.Rmd @@ -0,0 +1,263 @@ +--- +title: "Antimalarial Resistance" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Antimalarial Resistance} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + + +```{r setup} + +# Load the requisite packages: +library(malariasimulation) + +# Set colour palette: +cols <- c("#E69F00", "#56B4E9", "#009E73", "#CC79A7","#F0E442", "#0072B2", "#D55E00") + +``` + +## Introduction +One of the major threats to the continued success of efforts to reduce the burden of malaria is the evolution and spread of resistance to the antimalarial drugs used to treat uncomplicated cases of malaria. The most effective frontline antimalarials are a class of drugs called artemisinin combination therapies (ACTs). ACTs combine a fast-acting, short-lived artemisinin derivative (e.g. artemether), with a slower-acting, longer-lasting partner drug (e.g. lumefantrine) that clears the remaining parasites from the patient's system. Efforts to understand the effect of resistance to ACTs on malaria morbidity and mortality, and to develop strategies to control or mitigate the spread of resistance, would benefit from insights derived from mathematical modelling. Building on the model developed by Slater et al. (2016), `malariasimulation` provides the functionality to simulate the effects of resistance to the artemisinin and/or partner drug components to multiple, independent ACTs, on malaria transmission dynamics. + +Resistance to the artemisinin component of an ACT can result either in slow parasite clearance (SPC), in which treatment with an ACT takes longer than 3 days to fully clear patients with resistant parasites, or early treatment failure (ETF), in which the ACT fails to clear the infection and the individual develops a clinical infection. Resistance to the partner drug, where the partner drug fails to clear the parasite after the artemisinin derivative is depleted, results in infections recrudescing to either clinical (D) or asymptomatic infections (A). Resistance to the partner drug can also result in individuals developing a novel, resistant infection following treatment, as the prophylaxis provided by the ACT fails to protect the individual against reinfection by a resistant strain. In the following vignette, we illustrate how to parameterise and run `malariasimulation` simulations with resistance to ACTs deployed as a clinical treatment + +## Using set_antimalarial_resistance() to parameterise resistance +Simulations capturing the effects of resistance to clinical treatment using antimalarial drugs are parameterised using the `set_antimalarial_resistance()` function. This function appends user-defined resistance parameters to a `malariasimulation` parameter list and accepts ten inputs. The first is a list of `malariasimulation` parameters to append the resistance parameters to, and the second the index of the `drug` for which resistance is being parameterised, as set using the `set_drugs()` function. The `set_antimalarial_resistance()` function requires the `timesteps`, `artemisinin_resistance`, `partner_drug_resistance`, `slow_parasite_clearance_prob`, `early_treatment_failure_prob`, `late_clinical_failure_prob`, `late_parasitological_failure_prob`, and `reinfection_prob` inputs to be of equal length so that, for each time step in which an update occurs, a value is available for each parameter. Finally, the `slow_parasite_clearance_time` parameter represents the mean residence time, in days, for artemisinin-resistant individuals experiencing slow parasite clearance (SPC) in the Treated compartment, and must be input as a single, positive value. + +## Simulating static resistance +To illustrate how to parameterise resistance to an ACT using the `set_antimalarial_resistance()` function, we'll set-up and run three simulations. The first simulates malaria transmission in the absence of interventions or resistance. The second simulates a simple regime of clinical treatment in which 80% of clinical cases are treated with artemether lumefantrine (AL), beginning after one year, in the absence of antimalarial resistance. The third simulates the same clinical treatment programme but with resistance to the artemisinin component of AL emerging after two years. For illustrative purposes, we assume that the proportion of infections resistant to the artemisinin component of AL increases from 0% to 80%, and that these infections have a 90% chance of resulting in early treatment failure. + +### Parameterisation +First, we load the default `malariasimulation` model parameters, using the `overrides` argument to increase the human population. The human and mosquito population parameters are then calibrated to a user-specified initial EIR using the `set_equilibrium()` function. Next, we load the in-built parameters for the antimalarial drug AL and append them to the parameter list using `set_drugs()`. We can then use `set_clinical_treatment()` to specify a clinical treatment regime, beginning after one year, that treats, on average, 80% of the clinical cases of malaria with AL (`AL_params`). The `set_antimalarial_resistance()` function is then used to specify a scenario in which resistance is initially absent from the population, but after two years the proportion of malaria infections that are resistant to the artemisinin component of AL rises to 80%. We also set the probability that artemisinin-resistant infections result in early treatment failure to 0.9. In the current instance, we've set the proportion of infections resistant to the AL partner drug to 0% and the probabilities of other resistant infection outcomes to zero for simplicity. + +```{r, eval = TRUE} + +# Specify the time steps over which to simulate: +timesteps <- 365 * 3 + +# Specify an initial EIR to calibrate to: +init_EIR <- 8 + +# Specify a time step for treatment to begin: +treatment_start <- (1 * 365) + 1 + +# Specify a time step for resistance to emerge: +resistance_start <- (2 * 365) + 1 + +# Load the base malariasimulation parameter set: +simparams <- get_parameters( + overrides = list( + human_population = 10000) +) + +# Calibrate to the initial EIR: +simparams <- set_equilibrium(parameters = simparams, init_EIR = init_EIR) + +# Append the parameters for artemether lumefantrine (AL) to the parameter list: +simparams_clin_treatment <- set_drugs(parameters = simparams, drugs = list(AL_params)) + +# Use the set_clinical_treatment() function to specify a treatment regime for AL: +simparams_clin_treatment <- set_clinical_treatment(parameters = simparams_clin_treatment, + drug = 1, + timesteps = treatment_start, + coverages = 0.8) + +# Use the set_antimalarial_resistance() function to specify resistance to AL: +simparams_resistance <- set_antimalarial_resistance(parameters = simparams_clin_treatment, + drug = 1, + timesteps = c(0, resistance_start), + artemisinin_resistance = c(0, 0.8), + partner_drug_resistance = rep(0, 2), + slow_parasite_clearance_prob = rep(0, 2), + early_treatment_failure_prob = c(0, 0.9), + late_clinical_failure_prob = rep(0, 2), + late_parasitological_prob = rep(0, 2), + reinfection_prob = rep(0, 2), + slow_parasite_clearance_time = 10) + +``` + +### Simulation +We can now use the `run_simulation()` to simulate the three scenarios for which we have established parameter lists above. + +```{r, eval = TRUE} + +# Baseline: No-intervention, no resistance simulation: +sim_out_baseline <- run_simulation(timesteps = timesteps, parameters = simparams) + +# Clinical treatment with no antimalarial drug resistance: +sim_out_clin_treatment <- run_simulation(timesteps = timesteps, parameters = simparams_clin_treatment) + +# Clinical treatment with antimalarial drug resistance: +sim_out_resistance <- run_simulation(timesteps = timesteps, parameters = simparams_resistance) + +``` + +### Visualisation +With the outputs from the `run_simulation()` function, we can visualise the effect of resistance on malaria transmission by plotting the prevalence of malaria in children aged 2-10 years old (*Pf*PR~2-10~) over time. + +```{r, fig.align = 'center', eval = TRUE} + +# In each timestep, calculate PfPR_2-10 and add it to as a new column for each simulation output: +sim_out_baseline$pfpr210 = sim_out_baseline$n_detect_730_3650/sim_out_baseline$n_730_3650 +sim_out_clin_treatment$pfpr210 = sim_out_clin_treatment$n_detect_730_3650/sim_out_clin_treatment$n_730_3650 +sim_out_resistance$pfpr210 = sim_out_resistance$n_detect_730_3650/sim_out_resistance$n_730_3650 + +# Plot the prevalence through time in each of the three simulated scenarios: +plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE) +plot(x = sim_out_baseline$timestep, y = sim_out_baseline$pfpr210, + xlab = "Time (days)", + ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.7, + ylim = c(0, 1), type = "l", lwd = 2, xaxs = "i", yaxs = "i", + col = cols[3]) + +# Add a line for the clinical treatment scenario: +lines(x = sim_out_clin_treatment$timestep, + y = sim_out_clin_treatment$pfpr210, + col = cols[4]) + +# Add a line for the clinical treatment with resistance scenario: +lines(x = sim_out_resistance$timestep, + y = sim_out_resistance$pfpr210, + col = cols[7]) + +# Add lines to indicate the initiation of treatment and resistance: +abline(v = treatment_start, lty = "dashed") +abline(v = resistance_start, lty = "dashed") + +# Annotate the added vlines: +text(x = treatment_start + 10, y = 0.9, labels = "Start of\nTreatment", adj = 0, cex = 0.7) +text(x = resistance_start + 10, y = 0.9, labels = "Start of\nResistance", adj = 0, cex = 0.7) + +# Add gridlines: +grid(lty = 2, col = "grey80", nx = NULL, ny = NULL, lwd = 0.5); box() + +# Add a legend: +legend(x = 20, y = 0.99, legend = c("Baseline", "Treatment", "Resistance"), + col= c(cols[3:4], cols[7]), box.col = "white", + lwd = 1, lty = c(1, 1), cex = 0.7) + +``` + +In the absence of clinical treatment or resistance, prevalence is comparable between all three scenarios for the first year. Following the initiation of clinical treatment at the beginning of the second year, *Pf*PR~2-10~ approximately halves relative to the no-intervention baseline. However, following the introduction of artemisinin resistance at the beginning of the third year, early treatment failure causes the *Pf*PR~2-10~ to increase to an intermediate level in the resistance scenario. + +## Simulating dynamic resistance +We can also capture scenarios in which resistance to a drug changes through time. To illustrate, we'll establish and simulate a scenario in which resistance to sulfadoxine-pyrimethamine amodiaquine (SP-AQ) is absent from the population in the first year, but emerges in the third year and doubles in proportion each year thereafter until 100% of infections are artemisinin resistant. For simplicity, we'll assume only artemisinin resistance is present in the population, and resistance to artemisinin results only, and always, in early treatment failure. + +### Parameterisation +First, we store in vectors the artemisinin resistance proportions and the time steps on which they will be updated in the simulation. We also create a vector of early treatment failure probabilities which, for simplicity, we assume remain at 1 for each update. Next, we load the default `malariasimulation` parameter set, specifying a larger population size and seasonal transmission, and append the parameters for SP-AQ (`SP_AQ_params`) using the `set_drugs()` function. We'll specify a simple treatment regimen using `set_clinical_treatment()` where, on average, 80% of clinical cases are treated with SP-AQ, beginning after one year. We then specify a resistance schedule in which artemisinin resistance is introduced at a proportion of 0.2 after 3 years, and doubles each year thereafter until all infections are artemisinin resistant. Finally, we calibrate the human and mosquito population parameters to a defined entomological inoculation rate (EIR) and are ready to run the simulation. + +```{r, eval = TRUE} + +# Specify the time steps over which to simulate: +timesteps <- 365 * 8 + +# Set an initial EIR value: +initial_eir <- 8 + +# Specify the proportion of infections that are artemisinin resistant: +resistance_updates <- c(0, 0.2, 0.4, 0.8, 1) + +# Specify the time steps in which to update the resistance parameters: +resistance_update_timesteps <- c(0, seq(3*365, 6*365, by = 365)) + +# Specify the probability artemisinin resistance infections result in early treatment failure: +early_treatment_failure_updates <- rep(1, length(resistance_update_timesteps)) + +# Load the base malariasimulation parameter set, with seasonal transmission: +simparams <- get_parameters( + overrides = list( + human_population = 1000, + model_seasonality = TRUE, + g0 = 0.284596, + g = c(-0.317878,-0.0017527,0.116455), + h = c(-0.331361,0.293128,-0.0617547) +)) + +# Append the parameters for sulfadomamine pyrimethaine (SP-AQ) to the parameter list: +simparams <- set_drugs(parameters = simparams, drugs = list(SP_AQ_params)) + +# Use the set_clinical_treatment() function to specify a treatment regime for SP-AQ: +simparams <- set_clinical_treatment(parameters = simparams, + drug = 1, + timesteps = 365 * 1, + coverages = 0.8) + +# Use the set_antimalarial_resistance() function to specify resistance to SP-AQ: +simparams <- set_antimalarial_resistance(parameters = simparams, + drug = 1, + timesteps = resistance_update_timesteps, + artemisinin_resistance = resistance_updates, + partner_drug_resistance = rep(0, length(resistance_update_timesteps)), + slow_parasite_clearance_prob = rep(0, length(resistance_update_timesteps)), + early_treatment_failure_prob = early_treatment_failure_updates, + late_clinical_failure_prob = rep(0, length(resistance_update_timesteps)), + late_parasitological_prob = rep(0, length(resistance_update_timesteps)), + reinfection_prob = rep(0, length(resistance_update_timesteps)), + slow_parasite_clearance_time = 10) + +# Calibrate the parameters to an initial EIR: +simparams <- set_equilibrium(parameters = simparams, init_EIR = initial_eir) + +``` + +### Simulation +We can now use our parameter list to run the simulation using the the `run_simulation()` function. + +```{r, eval = TRUE} + +# Run the simulation: +dynamic_resistance_output <- run_simulation(timesteps = timesteps, parameters = simparams) + +``` + +### Visualisation +We can visualise the effect of increasing resistance through time by plotting the *Pf*PR~2-10~. We've added vertical lines to indicate when clinical treatment begins, and when the proportion of infections resistant to artemisinin is updated. + +```{r, fig.align = 'center', eval = TRUE} + +# Calculate the prevalence (PfPR210) through time: +dynamic_resistance_output$pfpr210 <- dynamic_resistance_output$n_detect_730_3650/dynamic_resistance_output$n_730_3650 + +# Open a new plotting window and add a grid: +plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE) +plot(x = dynamic_resistance_output$timestep, + y = dynamic_resistance_output$pfpr210, + xaxt = "n", + xlab = "Time (years)", + ylab = expression(paste(italic(Pf),"PR"[2-10])), cex = 0.8, + ylim = c(0, 1), type = "l", lwd = 2, xaxs = "i", yaxs = "i", + col = cols[3]) + +# Specify x-axis ticks and labels: +axis(1, at = seq(0, 8 * 365, by = 365), labels = seq(0, 8 * 365, by = 365)/365) + +# Add a line indicating the start of the clinical treatment: +abline(v = 365, lty = "dotted") + +# Add lines indicating when resistance is updated: +abline(v = resistance_update_timesteps, lty = "dashed") + +# Add a line highlighting the maximum PfPR_2-10 value prior to treatment or resistance: +abline(h = max(dynamic_resistance_output$pfpr210[1:365]), col = "red") + +# Add annotations for the vlines: +text(x = 365 + 30, y = 0.6, labels = "Clin. treat. begins", adj = 0, cex = 0.8, srt = 90) +text(x = resistance_update_timesteps[2:5] + 30, y = 0.6, labels = paste0("Art. Res. = ", resistance_updates[2:5]), + adj = 0, cex = 0.8, srt = 90) + +``` + +Looking at the figure, we can see that the *Pf*PR~2-10~ decreases over the two years following the onset of clinical treatment in the absence of artemisinin resistance. However, as resistance is introduced and increases through time, the *Pf*PR~2-10~ increases towards the pre-intervention seasonal peak as SP-AQ becomes increasingly ineffective in the treatment of clinical cases of malaria. + +## References +Slater, H.C., Griffin, J.T., Ghani, A.C. and Okell, L.C., 2016. Assessing the potential impact of artemisinin and partner drug resistance in sub-Saharan Africa. Malaria journal, 15(1), pp.1-11. \ No newline at end of file