changes needed for paper revision, added combination of interventions

SwissTPH · Jun 10, 2024 · fed7a5b · fed7a5b
1 parent 244c7b3
commit fed7a5b
Show file tree

Hide file tree

Showing 35 changed files with 1,266 additions and 809 deletions.
diff --git a/.gitignore b/.gitignore
@@ -10,4 +10,5 @@ bbl
 /Meta/
 .Rproj.user/
 .Rproj.user
-/.Rproj.user
+/.Rproj.user
+docs
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -30,7 +30,10 @@ Imports:
 Suggests: 
     knitr,
     rmarkdown,
-    testthat
+    testthat,
+    DiagrammeR, 
+    ggpubr,
+    Hmisc
 VignetteBuilder: knitr
 License: GPL (>=2)
 Encoding: UTF-8

diff --git a/NAMESPACE b/NAMESPACE
@@ -2,6 +2,8 @@
 
 export(AnophelesModel)
 export(build_model_obj)
+export(calculate_combined_impact)
+export(calculate_combined_impact_var)
 export(calculate_impact)
 export(calculate_impact_ip)
 export(calculate_impact_var)

diff --git a/R/AnophelesModel.r b/R/AnophelesModel.r
diff --git a/R/AnophelesModel_init.R b/R/AnophelesModel_init.R
@@ -116,6 +116,10 @@ def_vector_params = function(mosquito_species = "Anopheles gambiae",
 #' match one of the Anopheles species provided in the package.
 #' To see all available Anopheles species, use the function
 #' \code{list_all_species()}.
+#' @param vec_params list object created with the \code{def_vector_params()}
+#' function. The mosquito species attribute of this object needs to match the
+#' previously defined \code{mosquito_species} argument to this function.
+#' Default value is equal to the default \code{def_vector_params()}.
 #' @param host_table data frame with custom host-specific values.
 #' Must be provided if custom parameter values should be used instead of
 #' the ones in the package database.
@@ -158,7 +162,17 @@ def_vector_params = function(mosquito_species = "Anopheles gambiae",
 #' @export
 #'
 def_host_params = function(mosquito_species = "Anopheles gambiae",
+                           vec_params = def_vector_params(),
                            host_table = NULL, verbose = TRUE) {
+    # Check if the mosquito species is the same as the one in the vec_params
+    # First check if the given mosquito species is in the package database
+    if(mosquito_species != vec_params$species_name) {
+        err_msg = paste("Mosquito species name" , mosquito_species,
+                         "does not match with the species name from vec_params:",
+                         vec_params$species_name)
+        stop(err_msg)
+    }
+
     # Parameters specified by the user
     if(!is.null(host_table)) {
         # Check if input is correct
@@ -167,7 +181,7 @@ def_host_params = function(mosquito_species = "Anopheles gambiae",
         checkmate::assertDataFrame(host_table, any.missing = FALSE,
                                    col.names = "named", nrows = 2,
                                    add = assertCol)
-        # with the same column names as vec_ent_params
+        # with the same column names as host_ent_param
         checkmate::assertSubset(colnames(host_table), colnames(host_ent_param),
                                 empty.ok = FALSE, fmatch = FALSE)
         # and numeric, positive entries less than 1 except for the name of
@@ -202,6 +216,9 @@ def_host_params = function(mosquito_species = "Anopheles gambiae",
         host_param = as.list(host_p)
     }
 
+    # Check if the host parameters are well defined
+    hosts_p = calc_hosts_ent_params(vec_params, host_param, maxpop = 1000)
+
     return(host_param)
 }
 
@@ -334,6 +351,14 @@ def_activity_patterns = function(activity = "default_Anopheles_gambiae") {
 #' @export
 #'
 build_model_obj = function(vec_p, hosts_p, activity, total_pop) {
+    # Check if the mosquito species is the same in the vec_p and hosts_p
+    if(vec_p$species_name != unique(hosts_p$species_name)) {
+        err_msg = paste("Mosquito species name", vec_p$species_name,
+                        "from argument vec_p does not match with the species name from hosts_p",
+                        unique(hosts_p$species_name))
+        stop(err_msg)
+    }
+
     # Calculate the mosquito death rate for biting each host type and the host
     # availability rate
     assertNumeric(total_pop, lower = 1)
@@ -455,9 +480,7 @@ build_model_obj = function(vec_p, hosts_p, activity, total_pop) {
 #' LLIN insecticide type
 #' \item \code{LLIN_country:} only for LLINs, a string containing the country
 #' where the LLIN characteristics were measured
-#' \item \code{effects:} List of intervention survival and effects on the
-#' mosquito oviposition cycle (same format as defined above for the input
-#' attributes)
+#' \item \code{model_p:} model object provided as input
 #' }
 #'
 #' @author Monica Golumbeanu, \email{[email protected]}

diff --git a/R/internal_aux_functions.R b/R/internal_aux_functions.R
@@ -437,3 +437,16 @@ prepare_GVI_snippet = function(species, best_fit, param_name, id) {
     return(GVI_result)
 }
 
+# functions used for the combination of interventions
+thin_and_truncate = function(vec_in, thin, nips, redeploy = 1) {
+    vec_out = vec_in[c(TRUE, rep(FALSE, times = (thin - 1)))]
+    vec_out = rep(vec_out, redeploy)
+    return(vec_out[seq(1:nips)])
+}
+
+replace_effect = function(col1, df0, thin, nips) {
+    col2 = thin_and_truncate(t(df0[,2]), thin = thin, nips = nips)
+    col3 = thin_and_truncate(t(df0[,3]), thin = thin, nips = nips)
+    df1 = unname(cbind(col1, col2, col3))
+    return(df1)
+}
diff --git a/R/internal_entomological_model.R b/R/internal_entomological_model.R
@@ -29,6 +29,7 @@ calc_hosts_ent_params = function(vec_p, h_p, maxpop) {
     m = 1; n = 2
     zetai = c(rep(NA, m),rep(vec_p$zeta.3,n-m))
     PA = 1 - vec_p$A0
+
     # Assuming that there are only unprotected humans and non-human hosts
     # in the absence of interventions (eq. 14&15 Chitnis et al 2010) :
     PA1 = vec_p$A0 * vec_p$M * vec_p$Chi *
@@ -57,6 +58,11 @@ calc_hosts_ent_params = function(vec_p, h_p, maxpop) {
     # in the absence of interventions (eq. 27 Chitnis et al 2010)
     muvA = ((1 - (PA + sum(PAi[m:n]))) / (1 - PA)) * (-log(PA) / vec_p$td)
 
+    if(muvA < 0) {
+        stop("Simulated mosquito dynamics during the feeding cycle is not reflected by the mosquito bionomics parameters!
+        Consider increasing the values for PBi, PCi, PDi, PEi in the def_host_param() function call.",
+             call. = FALSE)
+    }
     h_p$muvA = muvA
     h_p$alphai = alphai
     return(h_p)

diff --git a/R/internal_intervention_init.R b/R/internal_intervention_init.R
@@ -97,6 +97,7 @@ init_intervention_effects = function(h_params, ni) {
 #   Kvi: proportion of susceptible mosquitoes that become infected
 #           after biting any host of type i
 #   alphai: availability to mosquitoes for host of type i
+#   model_p: parameters of the entomological model
 calc_interv_effects_db = function(interv_obj, model_p, ip) {
 
     # Initialize intervention object effects values
@@ -106,13 +107,21 @@ calc_interv_effects_db = function(interv_obj, model_p, ip) {
     db_interventions = list_intervention_models()
     if (interv_obj$id != "No intervention") {
         if (interv_obj$id %in% db_interventions$Intervention) {
+            # Extract the reference and the validation for the model
+            summary_tab = interventions_param$interventions_summary
+            interv_reference = summary_tab[which(summary_tab$Parameterisation == interv_obj$parameterisation), "Reference"]
+            interv_validation = summary_tab[which(summary_tab$Parameterisation == interv_obj$parameterisation), "Validation_ref"]
+            message(paste("Selected intervention parameterisation:", interv_obj$parameterisation))
+            message(paste("Reference for the selected parameterisation:", interv_reference))
+            message(paste("Previous validation for the selected parameterisation:", interv_validation))
             f_name = paste("calc_", interv_obj$id, "_p", sep="")
             # update the parameters by calling the intervention function
             interv_obj = do.call(f_name,
                                       list(int_obj = interv_obj,
                                            vec_params = model_p$vec_params,
                                            activity_cycles = model_p$activity,
                                            nips = ip))
+            interv_obj$model_p = model_p
         } else {
             err_msg = paste0("No intervention model available for",
                              interv_obj$id, ". To check available interventions

diff --git a/R/internal_intervention_models.R b/R/internal_intervention_models.R
@@ -10,7 +10,6 @@
 # but differs from the original implementation (for the Albimanus paper) in not
 # making use of the vector_specific_parameter vector or Pii
 calc_IRS_p = function(int_obj, vec_params, int_host, activity_cycles, nips) {
-    print("IRS")
     IRS_models_params = interventions_param$IRS_params
     int_summary = interventions_param$interventions_summary
     parameter_set = IRS_models_params[IRS_models_params$Parameterisation ==
@@ -67,6 +66,7 @@ calc_IRS_p = function(int_obj, vec_params, int_host, activity_cycles, nips) {
         }
         int_obj$effects[[paste0(name, "_decay")]] = points[[name]]
     }
+    int_obj$duration = duration
     return(int_obj)
 }
 
@@ -276,6 +276,7 @@ calc_LLINs_p = function(int_obj, vec_params, activity_cycles, nips) {
     int_obj$effects$PBi_decay = ip_vecs$PBi_decay
     int_obj$effects$PCi_decay = ip_vecs$PCi_decay
     int_obj$effects$survival = net_demog$survival
+    int_obj$duration = duration
     return(int_obj)
 }
 
@@ -287,6 +288,11 @@ calc_House_screening_p = function(int_obj, vec_params, activity_cycles, nips) {
     int_obj$effects$alphai[,1] = int_obj$effects$alphai[,2] *
         (1 - adjustment_for_location("alphai", indoor_outdoor,
                                      "House_screening") * 0.59)
+    d_idx = which(interventions_param$interventions_summary$Parameterisation ==
+                      int_obj$parameterisation)
+    duration = as.double(interventions_param$interventions_summary[d_idx,
+                                                                   "Duration"])
+    int_obj$duration = duration
     return(int_obj)
 }
 

diff --git a/data/host_ent_param.RData b/data/host_ent_param.RData
diff --git a/data/interventions_param.RData b/data/interventions_param.RData
diff --git a/data/vec_ent_param.RData b/data/vec_ent_param.RData
diff --git a/docs/authors.html b/docs/authors.html
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml
@@ -1,7 +1,7 @@
-pandoc: 3.1.1
+pandoc: '0'
 pkgdown: 2.0.7
 pkgdown_sha: ~
 articles:
   AnophelesModel: AnophelesModel.html
-last_built: 2023-10-17T20:04Z
+last_built: 2024-06-10T21:54Z
 
diff --git a/extdata/00_create_base_KEN.R b/extdata/00_create_base_KEN.R
@@ -28,10 +28,10 @@ library(OMAddons)
 library(stringr)
 library(AnophelesModel)
 
-# Function that creates a list with all the elements which are specific for the 
+# Function that creates a list with all the elements which are specific for the
 # base xml (placeholders, interventions, etc.)
 create_baseList = function(country_name, sim_start, versionnum) {
-  
+
   ## Basic xml skeleton
   baseList = list(
     # Mandatory
@@ -50,7 +50,7 @@ create_baseList = function(country_name, sim_start, versionnum) {
     # diagnostics = list(),
     model = list()
   )
-  
+
   # Create demography
   baseList = defineDemography(baseList,
                               name = country_name,
@@ -76,7 +76,7 @@ create_baseList = function(country_name, sim_start, versionnum) {
                                              0.004809187,
                                              0.00332171),
                               upperbound = c(1, seq(5, 85, by = 5)))
-  
+
   ## Create monitoring snippet
   baseList[["monitoring"]] = list(
     name = "Yearly Surveys",
@@ -100,13 +100,13 @@ create_baseList = function(country_name, sim_start, versionnum) {
       )
     ),
     surveys = monitoringSurveyTimesGen(detectionLimit = 100, startDate = "1999-01-01", #sim_start
-                                       endDate = "2025-01-01", 
-                                       interval = list(days = c(5), months = c(1:12), years = c(1998:2025)), 
+                                       endDate = "2025-01-01",
+                                       interval = list(days = c(5), months = c(1:12), years = c(1998:2025)),
                                        simStart = sim_start),
     ## surveyAgeGroupsGen will write thirdDimension table to cache, important for postprocessing
     ageGroup = surveyAgeGroupsGen(lowerbound = 0, upperbounds = c(1, 2, 5, 10, 100))
   )
-  
+
   ## Entomology section: MANDATORY
   seasonalityParameters = list(`Anopheles gambiae` = list(annualEIR="15",
                                                           input="EIR",
@@ -127,32 +127,32 @@ create_baseList = function(country_name, sim_start, versionnum) {
   # Define vector species in the simulation
   vectorSpecies = c("Anopheles gambiae")
   mosquitoParameters = mosquitoParameterization(vectorSpecies)
-  
+
   baseList$entomology$vector=list()
-  baseList = defineEntomology(baseList, seasonalityParameters, 
+  baseList = defineEntomology(baseList, seasonalityParameters,
                               mosquitoParameters)
-  
+
   # Begin interventions for humans
-  interventionList = list(LLIN_interv = list(id="LLINs", 
+  interventionList = list(LLIN_interv = list(id="LLINs",
                                              description="test LLIN",
                                              parameterisation="LLINs01",
                                              LLIN_type="Default",
                                              LLIN_insecticide="Default",
                                              LLIN_country="Kenya"))
-  
+
   vectorControlParameters = vectorControlParameterization(vectorSpecies,
                                                           interventionList)
   # Define LLIN intervention
   baseList$interventions$human = list()
   baseList = define_vector_control(baseList, vectorControlParameters)
-  
+
   ## Deployment section
-  
+
   # ITN deployment:
-  
+
   # At the moment the workflow is not flexible enough to allow custom GVI values.
   # So we need to replace the GVI effect values directly in the base xml with:
-  
+
   # $deterrency_snippet$GVI_xml_snippet
     # <GVI>
     #   <anophelesParams mosquito="Anopheles gambiae" propActive="1"/>
@@ -161,46 +161,46 @@ create_baseList = function(country_name, sim_start, versionnum) {
     #   <preprandialKillingEffect value="0"/>
     #   <postprandialKillingEffect value="0"/>
     # </GVI>
-  
+
   # $preprandial_snippet$GVI_xml_snippet
   # <GVI name="GVI_LLINs" id="GVI_LLINs_1">
   #   <anophelesParams mosquito="Anopheles gambiae" propActive="1"/>
   #   <decay L="1.86705167298283" function="weibull" k="1.72179371530776"/>
   #   <deterrency value="0"/>
   #   <preprandialKillingEffect value="0.826069186265107"/>
   #   <postprandialKillingEffect value="0"/>
-  #   </GVI> 
-  
+  #   </GVI>
+
   # $postprandial_snippet$GVI_xml_snippet
   # <GVI name="GVI_LLINs" id="GVI_LLINs_1">
   #   <anophelesParams mosquito="Anopheles gambiae" propActive="1"/>
   #   <decay L="1.3772267189462" function="weibull" k="1.01424385951852"/>
   #   <deterrency value="0"/>
   #   <preprandialKillingEffect value="0"/>
   #   <postprandialKillingEffect value="0.548152586557366"/>
-  #   </GVI> 
-    
-  baseList = deploy_IT(baseList = baseList, component = "LLIN_interv", 
-                       effects=c("deterrency", "preprandialKillingEffect", 
+  #   </GVI>
+
+  baseList = deploy_IT(baseList = baseList, component = "LLIN_interv",
+                       effects=c("deterrency", "preprandialKillingEffect",
                                  "postprandialKillingEffect"),
-                       coverage = 0.6, dates = c("2000-01-01"))
-  
+                       coverage = 0.4, dates = c("2000-01-01", "2002-01-01"))
+
   # Importation: MANDATORY
   baseList = define_importedInfections_compat(baseList = baseList, 10, time = 0)
-  
+
   # Health system
   # Write health system: MANDATORY
-  baseList = define_health_system(baseList = baseList, 
+  baseList = define_health_system(baseList = baseList,
                                   pSeekOfficialCareUncomplicated1 = 0.5,
                                   pSeekOfficialCareUncomplicated2 = 0.5)
-  
+
   ## Specify seed: MANDATORY
-  baseList = write_end_compat(baseList = baseList, 
+  baseList = write_end_compat(baseList = baseList,
                               seed = "@seed@", modelname = "base")
-  
+
   return(baseList)
 }
 
 # For testing
-# a = create_baseList("Kenya", "1918-01-01", 44L) 
+# a = create_baseList("Kenya", "1918-01-01", 44L)