Merge pull request #4 from mrc-ide/scene2

Scene2

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: scene
 Type: Package
 Title: Future Scenario Builder for the Malariaverse
-Version: 0.1.2
+Version: 1.0.0
 Authors@R: c(
     person(
       given = "Pete",
@@ -15,7 +15,7 @@ Description: Functionality to help define and build future scenarios for malaria
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 Remotes:
   mrc-ide/netz
 Depends: 
@@ -27,7 +27,6 @@ Imports:
     patchwork,
     rlang (>= 0.4.11),
     stats,
-    tibble,
     tidyr,
     zoo
 Suggests: 

R/modify.R

@@ -26,19 +26,22 @@ expand_interventions <- function(interventions, max_year, group_var){
 #' @param target New value of change point
 #'
 #' @export
-set_change_point <- function(interventions, sites, var, year, target){
+set_change_point <- function(interventions, var, year, target, sites = NULL){
 
-  template <- sites
-  template$year <- year
-  template$target <- target
+  if(is.null(sites)){
+    template <- data.frame(year = year)
+  } else {
+    template <- sites
+    template$year <- year
+  }
 
-  index <- index_df(interventions, template[,c(names(sites), "year")])
+  index <- index_df(interventions, template)
 
   if(!all(is.na(interventions[index, var]))){
     stop(paste("Trying to overwrite existing value in interventions:", var))
   }
 
-  interventions[index, var] <- template$target
+  interventions[index, var] <- target
 
   return(interventions)
 }
@@ -68,7 +71,7 @@ ever_used <- function(interventions, var, group_var){
 #' @param group_var Site grouping
 #'
 #' @export
-last_used <- function(interventions, var, group_var){
+last_used <- function(interventions, var, group_var = NULL){
   interventions <- interventions |>
     dplyr::filter(!is.na(.data[[var]])) |>
     dplyr::group_by(dplyr::across(dplyr::all_of(group_var))) |>
@@ -88,7 +91,7 @@ last_used <- function(interventions, var, group_var){
 #' @param group_var Site grouping
 #'
 #' @export
-linear_interpolate <- function(interventions, vars, group_var){
+linear_interpolate <- function(interventions, vars, group_var = NULL){
   interventions <- interventions |>
     dplyr::group_by(dplyr::across(dplyr::all_of(group_var))) |>
     dplyr::mutate(dplyr::across(dplyr::all_of(vars), \(x) zoo::na.approx(x, na.rm = FALSE))) |>
@@ -101,12 +104,13 @@ linear_interpolate <- function(interventions, vars, group_var){
 #' @param interventions Site file interventions section
 #' @param group_var Site grouping
 #' @param not Character string naming Interventions not to be extrapolated
+#' @param dir Direction ("down" or "up")
 #'
 #' @export
-fill_extrapolate <- function(interventions, group_var, not = "itn_input_dist"){
+fill_extrapolate <- function(interventions, group_var = NULL, not = "itn_input_dist", dir = "down"){
   f_interventions <- interventions |>
     dplyr::group_by(dplyr::across(dplyr::all_of(group_var))) |>
-    tidyr::fill(-dplyr::all_of(not), .direction = "down") |>
+    tidyr::fill(-dplyr::all_of(not), .direction = dir) |>
     dplyr::ungroup()
   return(f_interventions)
 }

R/plot.R

@@ -9,8 +9,8 @@ plot_interventions_combined <- function(
     interventions,
     population,
     group_var,
-    include = c("itn_use", "itn_input_dist", "fitted_usage", "tx_cov", "irs_cov", "rtss_cov", "smc_cov", "pmc_cov"),
-    labels = c("ITN usage", "ITN model input", "ITN model usage", "Treatment", "IRS", "RTSS", "SMC", "PMC"),
+    include = c("itn_use", "itn_input_dist", "predicted_use", "tx_cov", "irs_cov", "rtss_cov", "r21_cov", "smc_cov"),
+    labels = c("ITN usage", "ITN model input", "ITN model usage", "Treatment", "IRS", "RTSS", "R21", "SMC"),
     text = c(11, 11),
     facet_rows = 4,
     linewidth = 1
@@ -50,15 +50,16 @@ plot_interventions_combined <- function(
 #' @param labels Labels for variables in include
 #' @param text_size Size of font
 #' @param facet_rows Number of rows for faceted plot
+#' @param linewidth Line width
 #'
 #' @return Intervention plot
 #' @export
 plot_interventions <- function(
     interventions,
     population,
     group_var,
-    include = c("itn_use", "itn_input_dist", "fitted_usage", "tx_cov", "irs_cov", "rtss_cov", "smc_cov", "pmc_cov", "lsm_cov"),
-    labels = c("ITN usage", "ITN model input", "ITN model usage", "Treatment", "IRS", "RTSS", "SMC", "PMC", "LSM"),
+    include = c("itn_use", "itn_input_dist", "predicted_use", "tx_cov", "irs_cov", "rtss_cov", "r21_cov", "smc_cov"),
+    labels = c("ITN usage", "ITN model input", "ITN model usage", "Treatment", "IRS", "RTSS", "R21", "SMC"),
     text_size = 8,
     facet_rows = 4,
     linewidth = 1
@@ -73,7 +74,7 @@ plot_interventions <- function(
 
   pd <- interventions
   if(nrow(interventions) > nrow(unique(interventions[, group_var]))){
-    pd <- scene:::aggregate_df(df = interventions, groups = group_var, weighted_mean_cols = include, w = "par")
+    pd <- aggregate_df(df = interventions, groups = group_var, weighted_mean_cols = include, w = "par")
   }
 
   pd <- pd |>

R/utils.R

@@ -10,7 +10,7 @@ index_df <- function(df, filter_row) {
   indexed_df <- df |>
     dplyr::mutate(index = 1:dplyr::n())
 
-  filter_row <- filter_row[ , names(filter_row) %in% names(df)]
+  filter_row <- filter_row[ , names(filter_row) %in% names(df), drop = FALSE]
 
   row_index <- filter_row |>
     dplyr::left_join(indexed_df, by = names(filter_row), multiple = "all") |>

R/vector_control.R

@@ -17,31 +17,25 @@ add_future_net_dist <- function(interventions, group_var, off_year_max = 0.2, cy
       data_year = max(.data$year[!is.na(.data$itn_input_dist)]),
       # Infer a large scale net distribution year
       cycle_peak = .data$year[which.max(.data$itn_use[.data$year < .data$data_year])],
-      # Ensure previously fitted distributions remain fixed
-      dl = ifelse(is.na(.data$itn_input_dist), 0, .data$itn_input_dist),
-      du = ifelse(is.na(.data$itn_input_dist), off_year_max, .data$itn_input_dist),
+      # Ensure future distribution are cyclical
+      du = ifelse(is.na(.data$itn_input_dist), off_year_max, 1),
       du = ifelse(is.na(.data$itn_input_dist) & ((.data$year - .data$cycle_peak) %% cycle_period == 0), 1, .data$du),
-      # Estimate remaining distributions
-      itn_input_dist = netz::fit_usage(
-        target_usage = .data$itn_use,
-        target_usage_timesteps = 1 + (.data$year - min(.data$year) + 0.5) * 365,
-        distribution_timesteps =  1 + (.data$year - min(.data$year)) * 365,
-        distribution_init = pmax(pmin(.data$itn_use * 0.75, .data$du), .data$dl),
-        distribution_lower = .data$dl,
-        distribution_upper = .data$du,
-        control = list(xtol_rel = 0.01)
-      )$par
-    ) |>
-    # Add the resulting model usage
-    dplyr::mutate(
-      fitted_usage = netz::population_usage(
-        distribution = .data$itn_input_dist,
-        distribution_timesteps = 1 + (.data$year - min(.data$year)) * 365,
-        timesteps = (1 + max(.data$year) - min(.data$year)) * 365,
-        half_life = 5 * 365)[1 + (.data$year - min(.data$year) + 0.5) * 365]
+      itn_input_dist = netz::usage_to_model_distribution(
+        .data$itn_use,
+        1 + (.data$year - min(.data$year) + 0.5) * 365,
+        1 + (.data$year - min(.data$year)) * 365,
+        distribution_upper  = .data$du,
+        mean_retention = .data$mean_retention[1]
+      ),
+      fitted_usage = netz::model_distribution_to_usage(
+        1 + (.data$year - min(.data$year) + 0.5) * 365,
+        .data$itn_input_dist,
+        1 + (.data$year - min(.data$year)) * 365,
+        mean_retention = .data$mean_retention[1]
+      )
     ) |>
     dplyr::ungroup() |>
-    dplyr::select(-c("dl", "du", "cycle_peak", "data_year"))
+    dplyr::select(-c("du", "cycle_peak", "data_year"))
   return(interventions)
 }
 
@@ -80,3 +74,5 @@ link_vector_control_parameters <- function(interventions){
     dplyr::left_join(irs_parameters, by = "irs_insecticide")
   return(interventions)
 }
+
+utils::globalVariables("closest")

README.Rmd

@@ -16,7 +16,7 @@ knitr::opts_chunk$set(
 # scene <img src="man/figures/Scene.png" align="right" width=30% height=30% />
 
 <!-- badges: start -->
-[![Project Status: WIP – Initial development is in progress, but there has not yet been a stable, usable release suitable for the public.](https://www.repostatus.org/badges/latest/wip.svg)](https://www.repostatus.org/#wip)
+[![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
 [![R-CMD-check](https://github.com/mrc-ide/scene/workflows/R-CMD-check/badge.svg)](https://github.com/mrc-ide/scene/actions)
 [![Coverage status](https://codecov.io/gh/mrc-ide/peeps/branch/main/graph/badge.svg)](https://codecov.io/github/mrc-ide/scene)
 <!-- badges: end -->

README.md

@@ -5,9 +5,9 @@
 
 <!-- badges: start -->
 
-[![Project Status: WIP – Initial development is in progress, but there
-has not yet been a stable, usable release suitable for the
-public.](https://www.repostatus.org/badges/latest/wip.svg)](https://www.repostatus.org/#wip)
+[![Project Status: Active – The project has reached a stable, usable
+state and is being actively
+developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active)
 [![R-CMD-check](https://github.com/mrc-ide/scene/workflows/R-CMD-check/badge.svg)](https://github.com/mrc-ide/scene/actions)
 [![Coverage
 status](https://codecov.io/gh/mrc-ide/peeps/branch/main/graph/badge.svg)](https://codecov.io/github/mrc-ide/scene)

data-raw/create_example_site.R

@@ -1,10 +1,11 @@
 itn_use <- c(0, 0.1, 0.2, 0.4, 0.4)
 net_input <- round(
-  netz::fit_usage(
-    target_usage = c(0, 0.1, 0.2, 0.4, 0.4),
-    target_usage_timesteps = 365 * 0:4 + 180,
+  netz::usage_to_model_distribution(
+    usage = c(0, 0.1, 0.2, 0.4, 0.4),
+    usage_timesteps = 365 * 0:4 + 180,
     distribution_timesteps = 1 + 365 * 0:4
-  )$par, 2
+  ),
+  2
 )
 
 example_site <- list(
@@ -19,23 +20,27 @@ example_site <- list(
       year = 1:5,
       itn_use = itn_use,
       itn_input_dist = net_input,
+      mean_retention = 1000,
       tx_cov = c(0, 0.3, 0.4, 0.45, 0.5),
       irs_cov = 0,
       rtss_cov = 0,
       smc_cov = 0,
-      pmc_cov = 0
+      pmc_cov = 0,
+      lsm_cov = 0
     ),
     data.frame(
       country = "Eg",
       site = "B",
       year = 1:5,
       itn_use = itn_use,
       itn_input_dist = net_input,
+      mean_retention = 1000,
       tx_cov = c(0, 0.3, 0.4, 0.45, 0.5),
       irs_cov = 0,
       rtss_cov = 0,
       smc_cov = c(0, 0, 0.8, 0.8, 0.8),
-      pmc_cov = 0
+      pmc_cov = 0,
+      lsm_cov = 0
     )
   ),
   population = rbind(