lm pcr prevalence explicit

DESCRIPTION

@@ -69,15 +69,13 @@ LazyData: true
 Remotes:
   mrc-ide/malariaEquilibrium,
   mrc-ide/individual
-Additional_repositories:
-  https://mrc-ide.r-universe.dev
 Imports:
     individual (>= 0.1.16),
     malariaEquilibrium (>= 1.0.1),
     Rcpp,
     statmod,
     MASS,
-    dqrng (>= 0.3.2.2),
+    dqrng (>= 0.4),
     sitmo,
     BH,
     R6,

R/biting_process.R

@@ -210,9 +210,7 @@ simulate_bites <- function(
 calculate_eir <- function(species, solvers, variables, parameters, timestep) {
   a <- human_blood_meal_rate(species, variables, parameters, timestep)
   infectious <- calculate_infectious(species, solvers, variables, parameters)
-  age <- get_age(variables$birth$get_values(), timestep)
-  psi <- unique_biting_rate(age, parameters)
-  infectious * a * mean(psi)
+  infectious * a
 }
 
 effective_biting_rates <- function(a, .pi, p_bitten) {

R/model.R

@@ -29,10 +29,10 @@
 #'  * 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
+#'  * n_detect_lm (or pcr): number of humans with an infection detectable by microscopy (or pcr) between an inclusive age range at this timestep. This
 #' defaults to n_detect_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
+#'  * p_detect_lm (or pcr): the sum of probabilities of detection by microscopy (or pcr) between an
 #' inclusive age range at this timestep. This
 #' defaults to p_detect_730_3650. Other age ranges can be set with
 #' prevalence_rendering_min_ages and prevalence_rendering_max_ages parameters.

R/render.R

@@ -4,8 +4,11 @@ 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 light 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
@@ -32,7 +35,8 @@ create_prevelance_renderer <- function(
 
     clinically_detected <- state$get_index_of(c('Tr', 'D'))
     detected <- clinically_detected$copy()$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]]
@@ -43,17 +47,22 @@ create_prevelance_renderer <- function(
         timestep
       ) 
       renderer$render(
-        paste0('n_detect_', lower, '_', upper),
+        paste0('n_detect_lm_', lower, '_', upper),
         in_age$copy()$and(detected)$size(),
         timestep
       )
       renderer$render(
-        paste0('p_detect_', lower, '_', upper),
+        paste0('p_detect_lm_', lower, '_', upper),
         in_age$copy()$and(clinically_detected)$size() + sum(
           prob[bitset_index(asymptomatic, in_age)]
         ),
         timestep
       )
+      renderer$render(
+        paste0('n_detect_pcr_', lower, '_', upper),
+        pcr_detected$copy()$and(in_age)$size(),
+        timestep
+      )
     }
   }
 }

src/Random.cpp

@@ -60,7 +60,7 @@ void Random::prop_sample_bucket(
 
     // all probabilities are the same
     if (heavy == n) {
-        for (auto i = 0; i < size; ++i) {
+        for (size_t i = 0; i < size; ++i) {
             *result = (*rng)((uint64_t)n);
             ++result;
         }
@@ -121,9 +121,9 @@ void Random::prop_sample_bucket(
     }
 
     // sample
-    for (auto i = 0; i < size; ++i) {
+    for (size_t i = 0; i < size; ++i) {
         size_t bucket = (*rng)((uint64_t)n);
-        double acceptance = dqrng::uniform01((*rng)());
+        double acceptance = rng->uniform01();
         *result = (acceptance < dividing_probs[bucket]) ? bucket :
             alternative_index[bucket];
         ++result;

tests/testthat/test-biology.R

@@ -16,10 +16,7 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR', {
     vector_models <- parameterise_mosquito_models(parameters, 1)
     solvers <- parameterise_solvers(vector_models, parameters)
     estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0)
-    age <- get_age(variables$birth$get_values(), 0)
-    psi <- unique_biting_rate(age, parameters)
-    omega <- mean(psi)
-    mean(estimated_eir) / omega / population * 365
+    mean(estimated_eir) / population * 365
   })
 
   expect_equal(
@@ -44,10 +41,7 @@ test_that('total_M and EIR functions are consistent with equilibrium EIR (with h
     vector_models <- parameterise_mosquito_models(parameters, 1)
     solvers <- parameterise_solvers(vector_models, parameters)
     estimated_eir <- calculate_eir(1, solvers, variables, parameters, 0)
-    age <- get_age(variables$birth$get_values(), 0)
-    psi <- unique_biting_rate(age, parameters)
-    omega <- mean(psi)
-    mean(estimated_eir) / omega / population * 365
+    mean(estimated_eir) / population * 365
   })
 
   expect_equal(

tests/testthat/test-render.R

@@ -34,19 +34,27 @@ test_that('that default rendering works', {
   mockery::expect_args(
     renderer$render_mock(),
     2,
-    'n_detect_730_3650',
+    'n_detect_lm_730_3650',
     2,
     timestep
   )
 
   mockery::expect_args(
     renderer$render_mock(),
     3,
-    'p_detect_730_3650',
+    'p_detect_lm_730_3650',
     1.5,
     timestep
   )
 
+  mockery::expect_args(
+    renderer$render_mock(),
+    4,
+    'n_detect_pcr_730_3650',
+    3,
+    timestep
+  )
+
 })
 
 test_that('that default rendering works when no one is in the age range', {

vignettes/Antimalarial_Resistance.Rmd

@@ -107,9 +107,9 @@ With the outputs from the `run_simulation()` function, we can visualise the effe
 ```{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
+sim_out_baseline$pfpr210 = sim_out_baseline$n_detect_lm_730_3650/sim_out_baseline$n_730_3650
+sim_out_clin_treatment$pfpr210 = sim_out_clin_treatment$n_detect_lm_730_3650/sim_out_clin_treatment$n_730_3650
+sim_out_resistance$pfpr210 = sim_out_resistance$n_detect_lm_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)
@@ -225,7 +225,7 @@ We can visualise the effect of increasing resistance through time by plotting th
 ```{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
+dynamic_resistance_output$pfpr210 <- dynamic_resistance_output$n_detect_lm_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)
@@ -365,27 +365,27 @@ plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)
 
 # Plot malaria prevalence in 2-10 years through time:
 plot(x = simulation_outputs$base$timestep,
-     y = simulation_outputs$base$n_detect_730_3650/simulation_outputs$base$n_730_3650,
+     y = simulation_outputs$base$n_detect_lm_730_3650/simulation_outputs$base$n_730_3650,
      xlab = "Time (days)",
      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])
 
 # Add the dynamics for no-intervention simulation
 lines(x = simulation_outputs$treatment$timestep,
-      y = simulation_outputs$treatment$n_detect_730_3650/simulation_outputs$treatment$n_730_3650,
+      y = simulation_outputs$treatment$n_detect_lm_730_3650/simulation_outputs$treatment$n_730_3650,
       col = cols[4])
 
 lines(x = simulation_outputs$etf$timestep,
-      y = simulation_outputs$etf$n_detect_730_3650/simulation_outputs$etf$n_730_3650,
+      y = simulation_outputs$etf$n_detect_lm_730_3650/simulation_outputs$etf$n_730_3650,
       col = cols[5])
 
 lines(x = simulation_outputs$spc$timestep,
-      y = simulation_outputs$spc$n_detect_730_3650/simulation_outputs$spc$n_730_3650,
+      y = simulation_outputs$spc$n_detect_lm_730_3650/simulation_outputs$spc$n_730_3650,
       col = cols[6])
 
 lines(x = simulation_outputs$etf_spc$timestep,
-      y = simulation_outputs$etf_spc$n_detect_730_3650/simulation_outputs$etf_spc$n_730_3650,
+      y = simulation_outputs$etf_spc$n_detect_lm_730_3650/simulation_outputs$etf_spc$n_730_3650,
       col = cols[7])
 
 # Add vlines to indicate when SP-AQ were administered:
@@ -407,4 +407,4 @@ legend(x = 3000, y = 0.99, legend = c("Baseline", "Treatment", "ETF-only", "SPC-
 ```
 
 ## 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.
+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.

vignettes/Carrying-capacity.Rmd

@@ -66,10 +66,10 @@ and can run the simulations and compare the outputs
 ```{r, fig.width = 7, fig.height = 4, out.width='100%'}
 set.seed(123)
 s <- run_simulation(timesteps = timesteps, parameters = p)
-s$pfpr <- s$n_detect_730_3650 / s$n_730_3650
+s$pfpr <- s$n_detect_lm_730_3650 / s$n_730_3650
 set.seed(123)
 s_seasonal <- run_simulation(timesteps = timesteps, parameters = p_seasonal)
-s_seasonal$pfpr <- s_seasonal$n_detect_730_3650 / s_seasonal$n_730_3650
+s_seasonal$pfpr <- s_seasonal$n_detect_lm_730_3650 / s_seasonal$n_730_3650
 
 par(mfrow = c(1, 2))
 plot(s$EIR_gamb ~ s$timestep, t = "l", ylim = c(0, 200), xlab = "Time", ylab = "EIR")
@@ -106,7 +106,7 @@ p_lsm <- p |>
   )
 set.seed(123)
 s_lsm <- run_simulation(timesteps = timesteps, parameters = p_lsm)
-s_lsm$pfpr <- s_lsm$n_detect_730_3650 / s_lsm$n_730_3650
+s_lsm$pfpr <- s_lsm$n_detect_lm_730_3650 / s_lsm$n_730_3650
 
 par(mfrow = c(1, 2))
 plot(s$EIR_gamb ~ s$timestep, t = "l", ylim = c(0, 150), xlab = "Time", ylab = "EIR")
@@ -144,7 +144,7 @@ p_stephensi <- p_stephensi |>
 
 set.seed(123)
 s_stephensi <- run_simulation(timesteps = timesteps, parameters = p_stephensi)
-s_stephensi$pfpr <- s_stephensi$n_detect_730_3650 / s_stephensi$n_730_3650
+s_stephensi$pfpr <- s_stephensi$n_detect_lm_730_3650 / s_stephensi$n_730_3650
 
 par(mfrow = c(1, 2))
 plot(s_stephensi$EIR_gamb ~ s_stephensi$timestep,
@@ -171,7 +171,7 @@ p_flexible <- p |>
 
 set.seed(123)
 s_flexible <- run_simulation(timesteps = timesteps, parameters = p_flexible)
-s_flexible$pfpr <- s_flexible$n_detect_730_3650 / s_flexible$n_730_3650
+s_flexible$pfpr <- s_flexible$n_detect_lm_730_3650 / s_flexible$n_730_3650
 
 par(mfrow = c(1, 2))
 plot(s$EIR_gamb ~ s$timestep, t = "l", ylim = c(0, 150), xlab = "Time", ylab = "EIR")

vignettes/EIRprevmatch.Rmd

@@ -117,7 +117,7 @@ simparams <- set_clinical_treatment(parameters = simparams,
 
 Having established a set of `malariasimulation` parameters, we are now ready to run simulations. In the following code chunk, we'll run the `run_simulation()` function across a range of initial EIR values to generate sufficient points to fit a curve matching *Pf*PR~2-10~ to the initial EIR. For each initial EIR, we first use the `set_equilibrium()` to update the model parameter list with the human and vector population parameter values required to achieve the specified EIR at equilibrium. This updated parameter list is then used to run the simulation.
 
-The `run_simulation()` outputs an EIR per time step, per species, across the entire human population. We first convert these to get the number of infectious bites experienced, on average, by each individual across the final year across all vector species. Next, the average *Pf*PR~2-10~ across the final year of the simulation is calculated by dividing the total number of individuals aged 2-10 by the number (`n_730_3650`) of detectable cases of malaria in individuals aged 2-10 (`n_detect_730_3650`) on each day and calculating the mean of these values. Finally, initial EIR, output EIR, and *Pf*PR~2-10~ are stored in a data frame.
+The `run_simulation()` outputs an EIR per time step, per species, across the entire human population. We first convert these to get the number of infectious bites experienced, on average, by each individual across the final year across all vector species. Next, the average *Pf*PR~2-10~ across the final year of the simulation is calculated by dividing the total number of individuals aged 2-10 by the number (`n_730_3650`) of detectable cases of malaria in individuals aged 2-10 (`n_detect_lm_730_3650`) on each day and calculating the mean of these values. Finally, initial EIR, output EIR, and *Pf*PR~2-10~ are stored in a data frame.
 
 ```{r}
 # Establish a vector of initial EIR values to simulate over and generate matching
@@ -158,7 +158,7 @@ malSim_prev <- lapply(
     mean(
       output[
         output$timestep %in% seq(4 * 365, 5 * 365),
-        'n_detect_730_3650'
+        'n_detect_lm_730_3650'
       ] / output[
         output$timestep %in% seq(4 * 365, 5 * 365),
         'n_730_3650'
@@ -276,7 +276,7 @@ library(cali)
 summary_mean_pfpr_2_10 <- function (x) {
 
  # Calculate the PfPR2-10:
- prev_2_10 <- mean(x$n_detect_730_3650/x$n_730_3650)
+ prev_2_10 <- mean(x$n_detect_lm_730_3650/x$n_730_3650)
 
  # Return the calculated PfPR2-10:
  return(prev_2_10)
@@ -325,8 +325,8 @@ malsim_sim <- run_simulation(timesteps = (simparams_malsim$timesteps),
                              parameters = simparams_malsim)
 
 # Extract the PfPR2-10 values for the cali and malsim simulation outputs:
-cali_pfpr2_10 <- cali_sim$n_detect_730_3650 / cali_sim$n_730_3650
-malsim_pfpr2_10 <- malsim_sim$n_detect_730_3650 / malsim_sim$n_730_3650
+cali_pfpr2_10 <- cali_sim$n_detect_lm_730_3650 / cali_sim$n_730_3650
+malsim_pfpr2_10 <- malsim_sim$n_detect_lm_730_3650 / malsim_sim$n_730_3650
 
 # Store the PfPR2-10 in each time step for the two methods:
 df <- data.frame(timestep = seq(1, length(cali_pfpr2_10)),

vignettes/MDA.Rmd

@@ -143,15 +143,15 @@ plot.new(); par(mar = c(4, 4, 1, 1), new = TRUE)
 
 # Plot malaria prevalence in 2-10 years through time:
 plot(x = mda_output$timestep,
-     y = mda_output$n_detect_730_3650/mda_output$n_730_3650,
+     y = mda_output$n_detect_lm_730_3650/mda_output$n_730_3650,
      xlab = "Time (days)",
      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])
 
 # Add the dynamics for no-intervention simulation
 lines(x = noint_output$timestep,
-     y = noint_output$n_detect_730_3650/noint_output$n_730_3650,
+     y = noint_output$n_detect_lm_730_3650/noint_output$n_730_3650,
      col = cols[4])
 
 # Add vlines to indicate when SP-AQ were administered:
@@ -282,15 +282,15 @@ legend("topleft",
 plot.new(); par(new = TRUE, mar = c(4, 4, 1, 1))
 
 # Plot malaria prevalence in 2-10 years through time:
-plot(x = smc_output$timestep, y = smc_output$n_detect_730_3650/smc_output$n_730_3650,
+plot(x = smc_output$timestep, y = smc_output$n_detect_lm_730_3650/smc_output$n_730_3650,
      xlab = "Time (days)", 
      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])
 
 # Add the dynamics for no-intervention simulation
 lines(x = noint_output$timestep,
-     y = noint_output$n_detect_730_3650/noint_output$n_730_3650,
+     y = noint_output$n_detect_lm_730_3650/noint_output$n_730_3650,
      col = cols[4])
 
 # Add lines to indicate SMC events:
@@ -305,4 +305,4 @@ legend("topleft",
        legend = c("SMC", "No-Int"),
        col= c(cols[3:4]), box.col = "white",
        lwd = 1, lty = c(1, 1), cex = 0.8)
-```
+```

vignettes/Metapopulation.Rmd

@@ -107,7 +107,7 @@ output3$mix <- 'perfectly-mixed'
 output <- rbind(output1, output2, output3)
 
 # get mean PfPR 2-10 by year
-output$prev2to10 = output$p_detect_730_3650 / output$n_730_3650
+output$prev2to10 = output$p_detect_lm_730_3650 / output$n_730_3650
 output$year = ceiling(output$timestep / 365)
 output$mix = factor(output$mix, levels = c('isolated', 'semi-mixed', 'perfectly-mixed'))
 output <- aggregate(prev2to10 ~ mix + EIR + year, data = output, FUN = mean)