reports and docs

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: macpan2
 Title: Fast and Flexible Compartmental Modelling
-Version: 1.14.0
+Version: 2.0.0
 Authors@R: c(
   person("Steve Walker", email="swalk@mcmaster.ca", role=c("cre", "aut")),
   person("Weiguang Guan", role="aut"),

R/formula_list_generators.R

@@ -352,8 +352,8 @@ MockChangeModel = function() {
 ##' @examples
 ##' sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
 ##' sir
-##' sir |> mp_euler()             |> mp_expand()
-##' sir |> mp_rk4()               |> mp_expand()
+##' sir |> mp_euler() |> mp_expand()
+##' sir |> mp_rk4() |> mp_expand()
 ##' sir |> mp_euler_multinomial() |> mp_expand()
 ##' 
 ##' @name state_updates

R/mp_tmb_model_spec.R

@@ -299,9 +299,10 @@ must_save_time_args = function(formulas) {
 #' being evaluated. For example, expressions that generate stochasticity should
 #' be listed in \code{sim_exprs} because TMB objective functions must be
 #' continuous.
-#' @param state_update (experimental) Optional character vector for how to 
-#' update the state variables when it is relevant. Options include `"euler"`, 
-#' `"rk4"`, and `"euler_multinomial"`.
+#' @param state_update Optional character vector for how to update the state 
+#' variables when it is relevant. Can be the name of any 
+#' \code{\link{state_updates}} function with the `mp_` prefix removed (e.g., 
+#' `"euler"`, `"rk4"`, and `"euler_multinomial"`).
 #' 
 #' @examples
 #' spec = mp_tmb_model_spec(

R/tmb_model_editors.R

@@ -199,13 +199,34 @@ mp_tmb_delete = function(model
 #' @param incidence_name Name of the incidence variable to be transformed.
 #' @param report_prob Value to use for the reporting probability; the
 #' proportion of cases that get reported.
-#' @param mean_delay Mean of the Gamma distribution of reporting delay times.
+#' @param mean_delay Mean of the Gamma distribution of reporting delay times
+#' measured in time steps. For example, if the mean delay time in days is 14
+#' and each time step in the model corresponds to one week, then `mean_delay`
+#' should be set to `2`.
 #' @param cv_delay Coefficient of variation of the Gamma distribution of
-#' reporting delay times.
+#' reporting delay times, with time measured in time steps (see more on this
+#' in the description of the `mean_delay` argument immediately above).
 #' @param reports_name Name of the new reports variable.
 #' @param report_prob_name Name of the variable containing `report_prob`.
 #' @param mean_delay_name Name of the variable containing `mean_delay`.
 #' @param cv_delay_name Name of the variable containing `cv_delay`.
+#' @examples
+#' sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
+#' sir_delays = mp_tmb_insert_reports(sir
+#'   , incidence_name = "infection"
+#'   , report_prob = 0.1
+#'   , mean_delay = 14
+#'   , cv_delay = 0.25
+#'   , reports_name = "reports"
+#' )
+#' sir_delays |> mp_print_during()
+#' (sir_delays 
+#'   |> mp_simulator(
+#'       time_steps = 50L
+#'     , outputs = c("infection", "reports"), 50L)
+#'   |> mp_trajectory()
+#' )
+#' 
 #' @export
 mp_tmb_insert_reports = function(model
   , incidence_name
@@ -231,7 +252,7 @@ mp_tmb_insert_reports = function(model
   kernel_length = qgamma(0.95, shape, scale = scale) |> ceiling() |> as.integer()
 
   expressions = c(
-      sprintf("%s ~ pgamma(1:%s, 1/(%s), %s * (%s^2))", map$dist, kernel_length + 1L, cv_delay_name, mean_delay_name, cv_delay_name)
+      sprintf("%s ~ pgamma(1:%s, 1/(%s^2), %s * (%s^2))", map$dist, kernel_length + 1L, cv_delay_name, mean_delay_name, cv_delay_name)
     , sprintf("%s ~ %s[1:%s] - %s[0:%s]", map$delta, map$dist, kernel_length, map$dist, kernel_length - 1L)
     , sprintf("%s ~ %s * %s / sum(%s)", map$kernel, report_prob_name, map$delta, map$delta)
     , sprintf("%s ~ convolution(%s, %s)", reports_name, incidence_name, map$kernel)

inst/starter_models/shiver/README.Rmd

@@ -338,9 +338,9 @@ Now we update the specification to include these defaults.
 spec = mp_tmb_update(spec
   , default = list(
       N = N
-    , V = V0
-    , I = I0
-    , H = H0
+    , V = 0 # V0
+    , I = 1 # I0
+    , H = 0 # H0
     , a = a
     , b = a
     , alpha = alpha
@@ -369,6 +369,7 @@ shiver_calibrator = mp_tmb_calibrator(
   # we also want to estimate initial E
   , par = c("beta_v","beta_s","E", "sigma", "gamma_h") 
   , outputs = states
+  , time = mp_sim_bounds(-15, 90, "steps")
 )
 # print to check
 shiver_calibrator
@@ -658,8 +659,8 @@ prior_distributions = list(
   , logit_report_prob = mp_normal(qlogis(0.1), sd_par)
   , logit_p = mp_normal(qlogis(1/4), 4)
   , log_E_I_ratio = mp_normal(0, sd_par)
-  , log_I = mp_normal(log(I0), sd_state)
-  , log_H = mp_normal(log(H0), sd_state)
+  , log_I = mp_normal(log(1), sd_state)
+  , log_H = mp_normal(log(1), sd_state)
   , log_R = mp_normal(log(1), sd_state)
 )
 
@@ -669,6 +670,9 @@ dd = rbind(reported_hospitalizations, reported_cases)
 # calibrate
 shiver_calibrator = mp_tmb_calibrator(
     spec = (multi_traj_spec 
+      |> mp_tmb_update(
+        #default = list(log_I = 0, log_H = 0, log_R = 0)
+      )
       |> mp_hazard()
     )
     # row bind both observed data
@@ -684,6 +688,7 @@ shiver_calibrator = mp_tmb_calibrator(
     # a flexible model of time variation.
   , tv = mp_rbf("rbf_beta", 5, sparse_tol = 1e-8)
   , outputs = c(states, "reported_incidence", "beta")
+  , time = mp_sim_bounds(-30, 90, "steps")
 )
 ```
 
@@ -713,7 +718,7 @@ Our prior for `sigma` is similar to the posterior, but `gamma_h` seems to have b
       , fill = "red"
       , alpha = 0.3
     )
-   + geom_point(data = dd, aes(time, value))
+   + geom_point(data = mutate(dd, time = time + 30L), aes(time, value))
    + facet_wrap(vars(matrix), scales = 'free')
    + theme_bw()
 )

inst/starter_models/shiver/README.md

@@ -41,7 +41,7 @@ Vaccines are typically subject to resource constraints and distribution
 strategies might prioritize vaccinations for specific subpopulations,
 such as immunocompromised people, to reduce bad outcomes. We model this
 with a flow of susceptibles entering the vaccination class. This flow
-could be a fixed rate, i.e. a constant proportion of the population
+could be a fixed rate, i.e. a constant proportion of the population
 receives a vaccine each time step, but instead we wish to capture a more
 realistic vaccination rate by allowing it to vary (see [Variable
 Vaccination Rate](#variable-vaccination-rate)).
@@ -138,7 +138,7 @@ We can implement vaccine constraints by adding more model complexity.
 Resource limitations create an upper bound on the number of vaccines
 that can be administered to susceptibles per time step. There is also
 the constraint that we can only vaccinate, at most, the current number
-of susceptibles i.e. the vaccination rate can be at most 1. These
+of susceptibles i.e. the vaccination rate can be at most 1. These
 constraints naturally lead us to consider a variable vaccination rate
 $\phi(S(t))$, instead of vaccinating a fixed proportion $\phi > 0$ per
 time step.

inst/starter_models/sir/README.md

@@ -150,12 +150,12 @@ print(spec |> mp_expand())
 #> ---------------------
 #> Default values:
 #> ---------------------
-#>  matrix row col value
-#>    beta           0.2
-#>   gamma           0.1
-#>       N         100.0
-#>       I           1.0
-#>       R           0.0
+#>  quantity value
+#>      beta   0.2
+#>     gamma   0.1
+#>         N 100.0
+#>         I   1.0
+#>         R   0.0
 #> 
 #> ---------------------
 #> Before the simulation loop (t = 0):
@@ -165,7 +165,7 @@ print(spec |> mp_expand())
 #> ---------------------
 #> At every iteration of the simulation loop (t = 1 to T):
 #> ---------------------
-#> 1: infection ~ S * (I * beta/N)
+#> 1: infection ~ S * (beta * I/N)
 #> 2: recovery ~ I * (gamma)
 #> 3: S ~ S - infection
 #> 4: I ~ I + infection - recovery