Implement suggest boxes and transmission type options for multi-wrapp…

…ers (and class)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: IPDMR
 Title: Support Code for the Introduction to Practical Disease Modelling Course
-Version: 0.4.3-1
-Date: 2024-11-07
+Version: 0.5.0-1
+Date: 2024-11-08
 Description: A collection of functions and classes that are used to illustrate
   fundamental concepts of disease modelling as part of a physical course series
   taught jointly by the University of Copenhagen and University of Sydney. The

NAMESPACE

@@ -34,6 +34,8 @@ importFrom(methods,new)
 importFrom(pbapply,pblapply)
 importFrom(rlang,.data)
 importFrom(rlang,set_names)
+importFrom(stats,dgamma)
+importFrom(stats,qgamma)
 importFrom(stats,quantile)
 importFrom(stats,rmultinom)
 importFrom(tibble,tibble)

R/BetweenGroupClass.R

@@ -117,7 +117,14 @@ BetweenGroupClass <- R6::R6Class(
 
       ## TODO: C++ function if private$.allcpp
 
-      trans_b <- colSums(private$.beta_matrix * sapply(private$.groups, \(x) x$I))
+      if(private$.trans_between == "frequency"){
+        multby <- sapply(private$.groups, \(x){ x$I / x$N })
+      }else if(private$.trans_between == "density"){
+        multby <- sapply(private$.groups, \(x) x$I)
+      }else{
+        stop("Unrecognised private$.trans_between")
+      }
+      trans_b <- colSums(private$.beta_matrix * multby)
       for(i in seq_along(private$.groups)){
         private$.groups[[i]]$trans_external <- trans_b[i]
         private$.groups[[i]]$update(d_time)
@@ -206,6 +213,7 @@ BetweenGroupClass <- R6::R6Class(
     .ngroups = numeric(),
     .groups = list(),
     .time = numeric(),
+    .trans_between = "density",
 
     .dummy=NULL
   ),
@@ -229,6 +237,13 @@ BetweenGroupClass <- R6::R6Class(
       private$.time
     },
 
+    #' @field transmission_between the between-group transmission type (frequency or density)
+    transmission_between = function(value){
+      if(missing(value)) return(private$.trans_between)
+      value <- match.arg(value, c("frequency","density"))
+      private$.trans_between <- value
+    },
+
     #' @field state a data frame of the current state of each group
     state = function(){
       lapply(private$.groups, function(x) x$state) |>

R/model_wrappers.R

@@ -31,9 +31,9 @@ NULL
 
 #' @rdname group_models
 #' @export
-sir_det <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmission_type="frequency", d_time=1, max_time=100){
+sir_det <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmission_type=c("frequency","density"), d_time=1, max_time=100){
 
-  qassert(transmission_type, "S1")
+  transmission_type <- match.arg(transmission_type)
 
   model <- make_group(update_type="deterministic", numE=0, numI=1, numR=1)
   #model <- WithinGroupModel$new(model_type="sir", update_type="deterministic", transmission_type=transmission_type, d_time=d_time)
@@ -120,9 +120,9 @@ seir_det <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0
 
 #' @rdname group_models
 #' @export
-sir_stoc <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmission_type="frequency", d_time=1, max_time=100, iterations=1L){
+sir_stoc <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmission_type=c("frequency","density"), d_time=1, max_time=100, iterations=1L){
 
-  qassert(transmission_type, "S1")
+  transmission_type <- match.arg(transmission_type)
   qassert(iterations, "X1")
 
   model <- make_group(update_type="stochastic", numE=0, numI=1, numR=1)
@@ -159,7 +159,7 @@ sir_stoc <- function(S=99, I=1, R=0, beta=0.25, gamma=0.2, delta=0.05, transmiss
 
 #' @rdname group_models
 #' @export
-seir_stoc <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_type="frequency", d_time=1, max_time=100, iterations=1L){
+seir_stoc <- function(S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_type=c("frequency","density"), d_time=1, max_time=100, iterations=1L){
 
   #qassert(transmission_type, "S1")
   transmission_type <- match.arg(transmission_type)

R/multi_wrappers.R

@@ -22,6 +22,8 @@
 #' @param vacc the vaccination rate parameter per unit time (must be positive)
 #' @param repl the replacement rate parameter per unit time (must be positive)
 #' @param cull the targeted culling rate parameter per unit time (must be positive)
+#' @param transmission_within the within-group transmission type (frequency or density)
+#' @param transmission_between the between-group transmission type (frequency or density)
 #' @param iterations the number of iterations to run (stochastic models only)
 #' @param d_time the desired time step (delta time)
 #' @param max_time the desired maximum time point (must be greater than the time step)
@@ -37,9 +39,11 @@ NULL
 
 #' @rdname multi_models
 #' @export
-multi_seir_det <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, d_time=1, max_time=100){
+multi_seir_det <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_within=c("frequency","density"), transmission_between=c("density","frequency"), d_time=1, max_time=100){
 
-  output <- multi_wrapper("deterministic", n_groups=n_groups, beta_matrix=beta_matrix, S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega, gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull, d_time=d_time, max_time=max_time, iterations=1L)
+  transmission_within <- match.arg(transmission_within)
+  transmission_between <- match.arg(transmission_between)
+  output <- multi_wrapper("deterministic", n_groups=n_groups, beta_matrix=beta_matrix, S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega, gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull, transmission_within=transmission_within, transmission_between=transmission_between, d_time=d_time, max_time=max_time, iterations=1L)
 
   class(output) <- c("ipdmr_dm", class(output))
   attr(output, "plot_caption") <- str_c("deterministic; ", n_groups, " groups")
@@ -49,9 +53,11 @@ multi_seir_det <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L,
 
 #' @rdname multi_models
 #' @export
-multi_seir_stoc <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, d_time=1, max_time=100, iterations=1L){
+multi_seir_stoc <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_within=c("frequency","density"), transmission_between=c("density","frequency"), d_time=1, max_time=100, iterations=1L){
 
-  output <- multi_wrapper("stochastic", n_groups=n_groups, beta_matrix=beta_matrix, S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega, gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull, d_time=d_time, max_time=max_time, iterations=iterations)
+  transmission_within <- match.arg(transmission_within)
+  transmission_between <- match.arg(transmission_between)
+  output <- multi_wrapper("stochastic", n_groups=n_groups, beta_matrix=beta_matrix, S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega, gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull, transmission_within=transmission_within, transmission_between=transmission_between, d_time=d_time, max_time=max_time, iterations=iterations)
 
   class(output) <- c("ipdmr_sm", class(output))
   attr(output, "iterations") <- iterations
@@ -61,20 +67,21 @@ multi_seir_stoc <- function(n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L,
 }
 
 
-multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, d_time=1, max_time=100L, iterations=1L){
+multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=0, numE=3L, beta=0.25, omega=0.2, gamma=0.2, delta=0.05, vacc=0, repl=0, cull=0, transmission_within=c("frequency","density"), transmission_between=c("density","frequency"), d_time=1, max_time=100L, iterations=1L){
 
   ## Check arguments:
+  transmission_within <- match.arg(transmission_within)
+  transmission_between <- match.arg(transmission_between)
   qassert(n_groups, "X1(0,)")
   assert_matrix(beta_matrix, "numeric", any.missing=FALSE, nrows=n_groups, ncols=n_groups)
   qassert(iterations, "X1(0,)")
   qassert(d_time,"N1(0,)")
   qassert(max_time,"N1(0,)")
 
-  stop("switch to using BetweenGroupClass and add numI/numR (also to model_wrappers) and transmission_type")
-
   model_pars <- list(
     S=S, E=E, I=I, R=R, numE=numE, beta=beta, omega=omega,
     gamma=gamma, delta=delta, vacc=vacc, repl=repl, cull=cull,
+    transmission_type=transmission_within,
     group_number=seq_len(n_groups)
   )
 
@@ -91,9 +98,7 @@ multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=
     group_split() |>
     set_names(str_c("Group_", seq_len(n_groups) |> format() |> str_replace_all(" ", "0"))) |>
     lapply(\(x){
-      #md <- WithinGroupModel$new("seir", update_type, "frequency", numE=x$numE, d_time=d_time, group_number=x[["group_number"]])
-      md <- make_group(update_type=update_type, numE=x$numE, group_name=x[["group_number"]], model_type="SEIR", implementation="C++")
-      md$transmission_type <- "frequency"
+      md <- make_group(update_type=update_type, numE=x$numE, numI=1, numR=1, group_name=x[["group_number"]])
       x$numE <- NULL
       x$group_number <- NULL
       for(nm in names(x)){
@@ -107,6 +112,7 @@ multi_wrapper <- function(update_type, n_groups, beta_matrix, S=99, E=0, I=1, R=
   ##  Set up the model:
   model <- BetweenGroupClass$new(models)
   model$beta_matrix <- beta_matrix
+  model$transmission_between <- transmission_between
   model$save()
 
   ## Iterate:

R/suggest_boxes.R

@@ -4,49 +4,59 @@
 #' This is a utility function that takes a desired mean and standard deviation for
 #' the distribution of transit times through a disease compartment, and suggests
 #' a suitable rate and number of sub-compartments that matches the provided mean
-#' and variance as closely as possible.  If a d_time argument is provided, then
-#' the illustrative plot of the Erlang (gamma) distribution will be overlaid
-#' with a negative binomial.
+#' and variance as closely as possible.
 #'
-#' @param mean
-#' @param sd
-#' @param d_time
+#' @param mean the desired mean waiting time
+#' @param sd the desired standard deviation in waiting time
 #'
 #' @returns the function returns a ggplot object illustrating the suggestions,
 #' and prints the suggested number of sub-compartments and rates to screen
 #'
 #' @examples
-#' suggest_boxes(mean=5, sd=1)
-#' suggest_boxes(mean=5, sd=3, d_time=0.1)
+#' suggest_boxes(mean=5, sd=1.4)
+#'
+#' @importFrom stats dgamma qgamma
 #'
 #' @export
-suggest_boxes <- function(mean, sd, d_time){
+suggest_boxes <- function(mean, sd){
 
-  if(missing(d_time)) d_time <- NA_real_
+  qassert(mean, "N1(0,)")
+  qassert(sd, "N1(0,)")
 
-  stop("IMPLEMENT ME")
+  rate <- round(1/mean, 2)
 
   variance <- sd^2
   shape <- mean^2/variance^2
 
-  ## TODO: calculate 99% CI from qgamma as from/to
-  curve(dgamma(x, shape, scale=mean/shape), from=max(0,mean-5*sd), to=mean+5*sd)
-
-  options <- c(floor(shape), ceiling(shape), ceiling(shape)+1)
-  shapes <- unique(options[options>0])
-  rates <- round(1/mean, 2)
-  scales <- shapes/mean
+  options <- seq(floor(shape), ceiling(shape), by=1L)
+  opt_shapes <- unique(options[options>0])
+  vars <- sqrt(mean^2/opt_shapes)
 
-  yields <- paste0(shapes, " subcompartments with overall rate ", rates, ", which yields:\n\tmean = ", round(1/rates, 2), "  (desired: ", round(mean,2), ")\n\tsd = ", round(shapes^0.5*scales, 2), "  (desired: ", round(sd,2), ")\n", sep="")
+  yields <- paste0(opt_shapes, " subcompartments with overall rate ", rate, ", which yields:\n\tmean = ", round(1/rate, 2), "  (desired: ", round(mean,2), ")\n\tsd = ", round(sqrt(vars), 2), "  (desired: ", round(sd,2), ")\n", sep="")
   yields
 
-  if(length(shapes)==1L){
+  if(length(opt_shapes)==1L){
     cat("The closest match is ", yields, sep="")
     #cat("The closest match is Erlang(shape=", shapes, ", scale=", scales, "), which yields:\n\tmean = ", round(shapes/rates, 2), "  (desired: ", round(mean,2), ")\n\tsd = ", round(shapes^0.5/scales, 2), "  (desired: ", round(sd,2), ")\n", sep="")
   }else{
     cat("Possible matches are:\n", paste0("\t", str_replace_all(yields, "\n\t","\n\t\t")))
   }
 
-  stop("Add mandatory N property to within-gp and show/add freq-based for beta matrix, and allow dens-based for multi wrapper")
+  seq_along(opt_shapes) |>
+    lapply(\(x){
+      shape <- opt_shapes[x]
+      scale <- mean/shape
+      range <- qgamma(c(0.001,1-0.001), shape, scale=scale)
+      tibble(
+        `#Boxes: `=as.character(shape), Time=seq(range[1],range[2],length.out=1e3),
+        `rate: `=as.character(rate)
+      ) |>
+        mutate(Density = dgamma(Time,shape,scale=scale))
+    }) |>
+    bind_rows() |>
+    ggplot(aes(x=Time, y=Density, col=`#Boxes: `, lty=`rate: `)) +
+    geom_line() +
+    theme_light() +
+    theme(legend.position="bottom")
 
 }