Merge pull request #50 from phac-nml-phrsd/48-pop-renorm-contact-mats

support specification of population size by age

.Rbuildignore

@@ -8,3 +8,4 @@
 ^docs$
 ^pkgdown$
 ^\.github$
+^vignettes/articles$

.Rprofile

@@ -0,0 +1 @@
+source("renv/activate.R")

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: EPACmodel
 Title: Use the Early Pandemic Age-structured Compartmental (EPAC) model
-Version: 1.2.2
+Version: 1.3.0
 Authors@R: c( 
     person("Irena", "Papst", , "irena.papst@phac-aspc.gc.ca", role = c("aut", "cre")),
     person("Michael", "WZ Li", , role = c("aut"))) 
@@ -27,7 +27,8 @@ Suggests:
     tibble,
     tidyr,
     DT,
-    forcats
+    forcats,
+    devtools
 Remotes: canmod/macpan2@98fb990
 Config/testthat/edition: 3
 VignetteBuilder: knitr

R/helpers_five-year-age-groups.R

@@ -7,24 +7,32 @@
 
 #' Make contact matrix from Mistry et al. data
 #'
-#' @param age.group.lower Numeric vector. Lower bounds for desired age groups
+#' @param age.group.lower (optional) Numeric vector. Lower bounds for aggregated age groups, if desired.
 #' @template param_setting.weight
+#' @param pop.new (optional) Numeric vector. Specifies new population sizes by age with which to renormalize contact matrices. Must match up with age groups given in `age.group.lower` (if provided) or original Mistry age groups (one year age groups up to 84, with a single age group for 85+).
 #'
 #' @return A list with:
 #'  - `p.mat`: the row-normalized matrix of contact probabilities
 #'  - `c.hat`: the vector of row sums from the contact rate matrix, giving the average contact rate per age group; used to set up the transmission vector to maintain the contact balance condition
 #' @export
 mk_contact_pars <- function(
-    age.group.lower,
-    setting.weight
+    age.group.lower = NULL,
+    setting.weight,
+    pop.new = NULL
 ){
-
-  # Pull population tables if we need to aggregate age groups
-  if(!is.null(age.group.lower)){
-    pop.orig = mk_pop_table() # original population
-    pop.new = mk_pop_table(age.group.lower)
+  # Arg checks
+  if (!is.null(age.group.lower) & !is.null(pop.new)) {
+    if (length(age.group.lower) != length(pop.new)) stop("Population sizes given in `pop.new` must be compatible with age groups requested using `age.group.lower`. Currently, the lengths of these two vectors does not match.")
   }
 
+  if (is.null(names(setting.weight))) stop("`setting.weight` vector must be named. Names should be 'school', 'work', 'household', and 'community'.")
+
+  if (any(names(setting.weight) != c("school", "work", "household", "community"))) stop("Names of `setting.weight` should be 'school', 'work', 'household', and 'community'.")
+
+  # Pull population data
+  pop.orig = mk_pop_table() # original population
+  if(!is.null(age.group.lower)) pop.orig.age = mk_pop_table(age.group.lower) # original pops aggregated by age
+
   # Load component matrices (default with age groups 0, 1, ..., 83, 84+)
   setting.list = c("school", "work", "household", "community")
   c.mat.list = lapply(setting.list, function(setting){
@@ -65,7 +73,7 @@ mk_contact_pars <- function(
       |> dplyr::group_by(age_susceptible, age_infectious)
       |> dplyr::summarise(value = sum(value), .groups = "drop")
       # Attach new aggregated population counts and row-normalise
-      |> dplyr::left_join(pop.new, by = dplyr::join_by(age_susceptible == age_group))
+      |> dplyr::left_join(pop.orig.age, by = dplyr::join_by(age_susceptible == age_group))
       |> dplyr::transmute(
         age_susceptible,
         age_infectious,
@@ -95,6 +103,14 @@ mk_contact_pars <- function(
     c.mat.list[[setting]] * setting.weight[setting]
   }))
 
+  # Population re-normalize, if desired
+  if(!is.null(pop.new)){
+    # get population data corresponding to current c.mat indices
+    if (!is.null(age.group.lower)) pop.orig <- pop.orig.age
+    # renormalize each row
+    c.mat <- c.mat * pop.orig$count / pop.new
+  }
+
   # Return row-normalized matrix, as well as original row sums (avg contacts
   # by age group)
   c.hat = rowSums(c.mat)

R/helpers_hosp.R

@@ -237,7 +237,8 @@ make_pvec <- function(values, mats){
     # recalculate contact parameters
     values$contact.pars <- mk_contact_pars(
         age.group.lower = seq(0, 80, by = 5),
-        setting.weight = values$setting.weight
+        setting.weight = values$setting.weight,
+        pop.new = values$pop
     )
 
     # allow user to pass flow directly and avoid recalculation from
@@ -275,7 +276,8 @@ make_pvec <- function(values, mats){
         # recalculate new contact parameters
         values$contact.pars.new <- mk_contact_pars(
             age.group.lower = seq(0, 80, by = 5),
-            setting.weight = values$setting.weight.new
+            setting.weight = values$setting.weight.new,
+            pop.new = values$pop
         )
 
         out <- c(

inst/models/hosp/run-after-simulator.R

@@ -1,6 +1,6 @@
 # set up model simulator to accept updated parameter values
 pf <- (data.frame()
-|> add_to_pf("state", values$state)
+  |> add_to_pf("state", values$state)
   |> add_to_pf("flow", flow)
   |> add_to_pf("transmission.", transmission)
   |> add_to_pf("contact.", contact)
@@ -36,7 +36,7 @@ if (scenario.name == "change-contacts") {
 
   # set up to accept updated scenario parameters
   pf <- (pf
-  |> add_to_pf("contact_changepoints", contact_changepoints_to_fill)
+    |> add_to_pf("contact_changepoints", contact_changepoints_to_fill)
     |> add_to_pf("contact_values", contact_values_to_fill)
     |> add_to_pf("transmission_values", transmission_values_to_fill)
   )

inst/models/hosp/run-before-simulator.R

@@ -2,7 +2,8 @@
 if(is.null(values$contact.pars)){
   values$contact.pars = mk_contact_pars(
     age.group.lower = seq(0, 80, by = 5),
-    setting.weight = values$setting.weight
+    setting.weight = values$setting.weight,
+    pop.new = values$pop
   )
 }
 
@@ -15,7 +16,8 @@ contact <- values$contact.pars$p.mat
 if(scenario.name == "change-contacts" & is.null(values$contact.pars.new)){
   values$contact.pars.new = mk_contact_pars(
     age.group.lower = seq(0, 80, by = 5),
-    setting.weight = values$setting.weight.new
+    setting.weight = values$setting.weight.new,
+    pop.new = values$pop
   )
 }