Condensation

Condensation is the process of summarizing the concatenated_state_vector so that it can be compared with available data. This comparison is done quantitatively through a negative log-likelihood function that measures how different the condensed concatenated_state_vector is from the data. The user can add any of the following variables to the condensed dataset.

1. Any state variable – already computed in state and stored in concatenated_state_vector
2. Element-wise sum of any set of state variables – already computed in sp but not stored
3. Any time-varying rate matrix element – already computed in ratemat and stored in concatenated_ratemat_nonzeros
4. Element-wise product of any two variables of type 1-3 – not computed
5. Element-wise sum of any variable of type 4 – not computed
6. Lag-n differences of any variables of type 1-5 – not computed
7. Convolutions of any variables of type 1-5 with a gamma-density kernel (see section below on computing convolutions) – not computed

The following are computed in classic MacPan.

1. All the S boxes
2. The sum over all the S boxes
3. The sum over all the E boxes
4. The sum over all the I boxes
5. The sum over all the H and H2 boxes
6. The sum over all the ICUs and ICUd boxes
7. The sum over all R boxes
8. The lag-1 difference of the total of all X boxes (with NA at the beginning) – called hosp
9. The sum over all X boxes
10. The lag-1 difference of the total of all D boxes (with NA at the beginning) – called death
11. The sum over all D boxes
12. All the foi columns
13. Incidence:
    • compute foi times S box for each epi category
    • compute row sums over the resulting columns
    • compute the convolution (see section below on computing convolutions) of all of these columns (including the row sums column) – called report

Computing Convolutions

The convolution is based on a gamma-density kernel.

The shape parameter of this gamma density depends on params['c_delay_cv']. \[ \text{shape} = \frac{1}{\text{c_delay_cv}^2} \] The scale parameter depends on params['c_delay_mean'] and the shape \[ \text{scale} = \text{(c_delay_mean)} \text{(c_delay_cv)}^2 \] Compute the vector \(\gamma\) using the TMB pgamma function, which takes a vector q and scalar shape and scale arguments. \[ \gamma = \text{pgamma}(q, \text{shape}, \text{scale}) \]

where \(q\) is a constant integer vector that gets created in the following way in R.

qmax = ceiling(qgamma(0.95, shape, scale))
q = 1:qmax

We will avoid computing qmax in C++ because it could introduce a discontinuity that could cause automatic differentiation to fail. This difference between refactored and classic MacPan only has the potential to cause numerical differences whenever params['c_delay_cv'] and params['c_delay_mean'] are varying – typically because they are in opt_pars – but I have not yet seen any examples where this is the case. To guarantee matched results at least for the initial parameter vector we will set qmax as follows on the R-side and pass it through a DATA_ macro.

qmax = ceiling(qgamma(0.95, 1/params[["c_delay_cv"]]^2, params[["c_delay_mean"]] * params[["c_delay_cv"]]^2))`

The delay kernel for the convolution is a numeric vector given by the following expression. \[ \kappa_i = (\text{c_prop})\frac{\delta_i}{\sum_j\delta_j} \] where \(\delta_j = \gamma_{j+1} - \gamma_j\) and \(\text{c_prop}\) is a scalar in the parameter vector.

The convolution, \(z\), of the time series of any state variable, \(x\), is given by the following. \[ z_{i+m-1} = \sum_{j = 1}^m \kappa_{m-j+1} x_{i+j-1} \] where \(i = 1,...,n-m\), \(m\) is the length of \(\kappa\), and \(n\) is the number of time steps (length of the input vector \(x\)). Note that the first \(m-1\) elements of \(z\) remain undefined and so they do not need to be stored on the C++-side because these \(m-1\) missing values can be added on the R-side.

Condensation Algorithm

In Progress: in collaboration with Weiguang

Condensation is the process of summarizing the history of the simulations. Before solving the condensation problem we should refactor how this history is currently summarized. Currently the concatenated_state_vector and concatenated_ratemat_nonzeros contain some history information, but it will be awkward to continue with this strategy. Instead we will introduce a matrix called simulation_history with rows corresponding to simulation steps and columns corresponding to the following (in the following order).

1. State variables
2. Changing rate matrix elements (there is one such element in this model)
3. Sums
4. Factrs
5. Element-wise expressions involving the first four items
6. Lag-n differences of the first five items
7. Convolutions of the first five items

Items 1-4 are computed at simulation time, and simply need to be added to the simulation_history matrix at each iteration.

Items 5-7 should be computed at the end of the simulation. These columns should be computed in the same order as they will appear in the simulation_history matrix. The input for these computations (5-7) will be stored in the following new data structures.

• Item 5 – element-wise expressions:
  • sr_count – vector the same length as the number of columns added in step 3a giving the number of simulation_history columns to use as input for each new column (note that input columns can be used more than once – see count vector for a similar structure)
  • sri – indices into the columns of simulation_history to use as input (see spi vector for a similar structure)
  • sr_modifier – bitwise integer encoding of the parse tree (see modifier for a similar data structure)
• Item 6 – lag-n differences:
  • lag_diff_sri – indices into the columns of simulation_history to use as input to lag-n difference operations
  • lag_diff_delay_n – vector same length as lag_diff_sri containing the delay (i.e. the “n” in lag-n) of the lag for each lag-n difference operation
• Item 7 – convolutions:
  • conv_sri – indices into the columns of simulation_history to use as input to convolution operations
  • conv_c_prop_idx, conv_c_delay_cv_idx, conv_c_mean_cv_idx – three vectors each the same length as conv_sri giving indices into params for extracting the parameters of the convolutions
  • conv_qmax – vector the same length as conv_sri giving the maximum integer in the q vector for each convolution

Add simulation_history to the report. On the R-side we could start using simulation_history where we currently use concatenated_* vectors, and eventually we could drop them.

Negative Binomial Negative Log Likelihood

We compare the negative binomial negative log likelihood with observed_vector as data and condensed_state_vector as predictions. Apologies that the word ‘negative’ is being used in two different ways – the first refers to the negative binomial distribution, and the second to the fact that we are multiplying the log likelihood by negative one. In R we would write this as follows.

-1 * dnbinom(observed_vector, mu = condensed_state_vector, size = nb_disp, log = TRUE)

In TMB, this negative binomial density function should be used for this (I think … how does size on the R-side correspond to var in TMB?).

The first three arguments to this function are vectors of the same size.

• observed_vector is an observed data stream that is provided by the user
• condensed_state_vector is described in detail above
• nb_disp is a vector of dispersion parameters, which is composed of values passed in through PARAMETER macros (see below for a description of how nb_disp is constructed)

The sum of these negative log likelihood values gives a loss function to be returned to the user.

TODO: how to handle clamping when we leave the support of the density? concerned about discontinuities.

Dispersion Parameters

We require an additional set of parameters to model negative binomial dispersion, nb_disp. The user has two options.

1. one dispersion parameter to be used in every element of nb_disp
2. one dispersion parameter is passed in for every variable in condensed_state_vector

Optional Prior Distribution Penalties

The user is allowed to specify a prior distribution for any element of params, tv_mult, or nb_disp. Initially only univariate normal prior distributions will be allowed, but this should be built to be extensible allowing for other distributions to be added in the future.

For each parameter with a prior, the user specifies the distribution (initially normal will be the only choice), and a parameter vector for the arguments of the corresponding density function. At each evaluation of the optimizer the log of the prior densities are subtracted from the negative log-likelihood function.

The TMB normal density function should be used. We should be extensible to allow future developers to add other distributions from here.

Optional Inverse Parameter Transformations

The user is allowed to express any value of params or tv_mult on one of several transformed scales. The corresponding inverse transformation is applied immediately as the parameter is read through the appropriate macro on the TMB-side. The user should be allowed to specify one of the following transformations.

• log10
• log
• logit

The corresponding inverse transformations that would be applied after the macros are the following.

• 10^x
• exp(x)
• 1/(1+exp(-x))

General Objective Function

Let \(l_{\psi_i}\left(x_i(\theta); y_i\right)\) be a loss function where \(x_i(\theta)\) is a simulated element of the history matrix, \(y_i\) is an associated observed data point, and \(\psi_i\) is a vector of loss function parameters. Note that we write the simulation values as functions of \(\theta\) to indicate their dependence on parameters. Every element \(i\) that belongs to the same column of the history matrix has the same value for \(\psi_i\), but to simplify notation \(\psi\) is sub-scripted by \(i\). Let \(r_{\phi_j}(\theta_j)\) be an optional regularization function for each parameter \(\theta_j\) being optimized, where \(\phi_j\) is a vector of regularization function parameters associated with parameter \(j\). The objective function is then given by the following.

\[ L_{\psi, \phi}(\theta; y) = \sum_i l_{\psi_i}\left(x_i(\theta); y_i\right) + \sum_j r_{\phi_j}(\theta_j) \] Currently the only objective function that is offered is the negative log negative binomial density, but the infrastructure we build assumes that other loss functions will be allowed in the future. Similarly the only regularizing function that is available is the negative log normal density.

Each parameter, \(\theta_j\), can be optionally passed into the objective function on a transformed scale. These transformations allow users to enforce limits on the domain of the parameter space (e.g. a logit transformation will keep parameter values between 0 and 1). The regularization functions are applied on the transformed scale, but the elements of the rate matrix (and therefore the simulations, \(x\)) are computed using the parameters on their original scale. Therefore, we have the following ordering of operations.

1. Parameter, \(\theta_j\), is passed in to the objective function on the transformed scale, \(f(\theta_j)\)
2. The regularization function, \(r_{\phi_j}(\theta_j)\), is computed and saved
3. The inverse transformation is applied to the parameter, \(\theta_j\), and the result is used to overwrite the previous value in the c++ params vector
4. The simulations are made
5. The objective function, \(L\), is computed, and this computation includes adding the saved values of the regularization functions
6. The objective function, \(L\), is returned by the c++ function objective_function

Available Loss Functions

• Negative log negative binomial density
  • Loss function parameters: dispersion
  • In TMB use dnbinom2 to compute the log negative binomial density
  • The arguments for dnbinom2 can be computed as follows:
    • x – observed data
    • mu – simulated data
    • var – mu + mu^2 / disp, where disp is the dispersion parameter
    • give_log – 1
  • Lower bound on the simulated values
    • Before passing the simulated data to the mu argument, the user may optionally choose to smoothly constrain the simulated data to be greater than a lower bound, eps, by applying the function mu_constrained = mu + eps*exp(-mu/eps)

Available Regularization Functions

• No regularization
  • Name to use in C++: flat
  • Regularization function parameters: initial value (\(\theta_\text{init}\))
  • \(r_{\theta_\text{init}}(\theta) = 0\)
  • The fact that the regularization parameter, \(\theta_\text{init}\), has no effect on the regularization function is an unusual special case. When regularization is applied, the ‘location’ of the regularization function is used as the initial value for the optimization. In the case of no regularization, we are using the regularization parameters to contain the initial value. One interpretation of this regularization parameter is as the ‘peak’ of a flat prior.
• Negative log of the normal density
  • Name to use in C++: normal
  • Regularization function parameters: mean (\(\bar{\theta}\)), standard deviation (\(\sigma_\theta\))
  • \(r_{\mu_f, \sigma_f}(\theta) = -\log\left[\mathcal{N}(f(\theta), \mu_f, \sigma_f)\right]\), where \(\mathcal{N}\) is the normal density
  • In TMB the dnorm function should be used with the following arguments
    • x – \(f(\theta)\) – value of the parameter on the transformed scale
    • mean – \(\mu_{f(\theta)}\) – first regularization parameter
    • sd – \(\sigma_{f(\theta)}\) – second regularization parameter
    • give_log ` – 1

The first regularization parameter for every regularization function is used as the starting value for the optimizer on the scale of the transformed parameter.

Available Parameter Transformations

• identity
  • index = 1
  • trans = \(f(x) = x\)
  • inverse = \(f^{-1}(x) = x\)
• log
  • index = 2
  • trans = \(f(x) = \log(x)\)
  • inverse = \(f^{-1}(x) = e^x\)
• log10
  • index = 3
  • trans = \(f(x) = log_{10}(x)\)
  • inverse = \(f^{-1}(x) = 10^x\)
• logit
  • index = 4
  • trans = \(f(x) = \log\left[\frac{x}{1-x}\right]\)
  • inverse = \(f^{-1}(x) = \frac{1}{1 + e^{-x}}\)
• cloglog
  • index = 5
  • trans = \(f(x) = \log\left[-\log(1-x)\right]\)
  • inverse = \(f^{-1}(x) = 1 - \exp\left[-\exp(x)\right]\)
• inverse
  • index = 6
  • trans = \(f(x) = \frac{1}{x}\)
  • inverse = \(f^{-1}(x) = \frac{1}{x}\)

General Objective Function

We will pass the following additional data through to TMB to specify the objective function.

• Variables: the following three integer vectors have one element per variable in the objective function and identify what loss function is used for each variable, as well as how many parameters those loss functions require
  • obs_var_id – id giving the column in the history matrix associated with each variable in the objective function
  • obs_loss_id – identifier of the loss function being used for each variable (currently only one – negative log negative binomial density)
  • obs_loss_param_count – number of loss function parameters associated with each variable
• Loss Function Parameters: the following vector has one element per loss function parameter (these are the \(\psi_i\) parameters described above in the general objective function section)
  • obs_spi_loss_param – index into sp containing the loss function parameters. the number of such parameters is sum(obs_loss_param_count)
• Observed Data the following three vectors have one element for each observed data point, and identify what elements in the history matrix should be compared with each observed value (each element of these vectors corresponds to each \(i\) in the description of the general objective function above)
  • obs_time_step – row in the history matrix
  • obs_history_col_id – column in the history matrix (match with obs_var_id above to determine the obs_loss_id and obs_loss_param_count)
  • obs_value – value observed in data to compare with the associated element in the history matrix
  • obs_do_sim_constraint – Boolean determining if the lower bound in obs_sim_lower_bound should be applied to the simulated data with smooth constraint function mu_constrained = mu + eps*exp(-mu/eps)
  • obs_sim_lower_bound – lower bound applied if obs_do_sim_constraint is true
• Optimized Parameters: the following vectors have one element per parameter that is being optimized
  • opt_param_id – indices into the params vector giving the parameters to be optimized
  • opt_trans_id – index identifying a transformation – see above for correspondence between indices and transformations
  • opt_count_reg_params – number of regularization parameters associated with each parameter to be optimized
  • opt_reg_family_id – index identifying a family of regularization functions – currently only two are available (flat, normal)
• Optimized Time-Varying Parameters (not yet in scope): similar to optimized parameters but for tv_mult
  • opt_tv_param_id – indices into the tv_mult vector giving the time varying multipliers to be optimized
  • opt_tv_trans_id – index identifying a transformation – see above for correspondence between indices and transformations
  • opt_tv_count_reg_params – number of regularization parameters associated with each parameter to be optimized
  • opt_tv_reg_family_id – index identifying a family of regularization functions – currently only two are available (flat, normal)
• Regularization Parameters:
  • opt_reg_params – regularization parameters – the size of this vector is opt_count_reg_params.sum()
  • opt_tv_reg_params – regularization parameters for tv_mult – the size of this vector is opt_tv_count_reg_params.sum()

Generalization of Outflow Vector

Let \(f_i^{(out)}\) be the \(i\)th element of the outflow vector, which is computed as follows. \[ f_i^{(out)} = \sum_{j \in \Omega_i}F_{ij} \]

Where \(\Omega_i\) is a subset of the indices into the state vector. Typically there are only a small number of unique index sets, with one \(\Omega_i\) for many values of \(i\) – for this reason the data structure for these indices only need to track each unique index set.

We can express the previous definition of the outflow vector (TODO: link to previous definition) in the terms of this generalization as a constant \(\Omega = \Omega_i, \forall i\).

Make State Vector

In previous spec versions, the C++ side needs to get the initial state vector directly from R. The difficulty with this strategy is that the initial state vector depends on the parameter vector, so when a new parameter vector gets updated in an optimization or ensemble run it becomes necessary to drop back to R causing a substantial performance bottleneck. So here we add the ability to create the initial state vector from the parameter vector in C++.

What follows is a detailed description of the procedure, which makes reference to new items in the data structure (see Data Structure below).

0. We first check to see if a state vector has been passed – if not, do the following steps to create a new state vector
1. Initialize two state vectors with all zeros
   • state – initializing state vector for the full non-linearized model
   • lin_state – initializing state vector for the linearized model
2. Make a copy (lin_params) of the parameter vector (params)
3. Replace some user-specified elements of lin_params with user supplied constant values (lin_param_vals, lin_param_count, lin_param_idx – these can be found in model$tmb_indices$linearized_params)
4. Replace some user-specified elements of lin_state with elements of lin_params (df_state_par_idx, df_state_count, df_state_idx – these can be found in model$tmb_indices$disease_free)
5. Compute the Jacobian matrix (jac), which requires the following steps in a function that takes lin_state and returns an updated lin_state
   • Create a rate matrix (lin_ratemat) from lin_state and lin_params using the existing make_ratemat_indices
   • Create a flow matrix (lin_flow) from lin_state and lin_ratemat
   • Create inflow vector (lin_inflow) from lin_flow by taking column sums
   • Create outflow vector (lin_outflow) from lin_flow by the method above – Generalization of Outflow Vector (linearized_outflow_row_count, linearized_outflow_col_count, linearized_outflow_rows, linearized_outflow_cols)
   • Update lin_state using lin_inflow and lin_outflow
6. Remove certain user-defined rows, columns, and elements from jac and lin_state to create a new vector eigvec (all_drop_eigen_idx – can be found in model$tmb_indices$initialization_mapping)
7. Compute the dominant eigenvector, eigvec of jac
8. Remove elements of eigvec to create eig_infected (eigen_drop_infected_idx – can be found in model$tmb_indices$initialization_mapping)
9. Normalize eig_infected to have elements that sum to one
10. Update state[all_to_infected_idx] = eig_infected * param[infected_idx] (infected_idx can be found in model$tmb_indices$initial_population)
11. Update state[susceptible_idx] = (1/susceptible_idx.size()) * (param[total_idx] - param[infected_idx]) (infected_idx and total_idx can be found in model$tmb_indices$initial_population)

Generalization of Outflow Vector

Objects of class flexmodel contains a new element called tmb_indices$outflow, with the following elements

• row_count –
• col_count –
• rows –
• cols –

These indices encode information in the \(\omega_i\) index sets (see Mathematical Description above).

This outflow item means that par_accum_indices is no longer required. Given that we do not need to worry about back compatibility yet, it will be dropped from this and all future spec versions.

Make State Vector

Objects of class flexmodel contains the following new elements.

• linearized_params – indices to update some elements of lin_param before computing the Jacobian of the linearized model
  • lin_param_vals – values used to fill the elements of lin_param
  • lin_param_count – the number of elements of lin_param that will be updated with each value in lin_param_vals
  • lin_param_idx – indices of lin_param to update
• disease_free – indices to initialize lin_state to be nearly disease-free, by setting some of its elements equal to some elements of lin_param
  • df_state_par_idx – indices of lin_param used to fill the elements of lin_state
  • df_state_count – the number of elements of lin_state that will be updated with each parameter identified by df_state_par_ix
  • df_state_idx – indices of lin_state to update
• linearized_outflow – same as outflow but for the linearized model
• initialization_mapping – indices into various subsets of the state vector and eigenvector that could be useful for initializing the state – see justification below this list
  • all_to_eigen_idx – indices into state and lin_state for identifying states in eigvec
  • all_to_infected_idx – indices into state and lin_state for identifying states in eig_infected
  • eigen_to_infected_idx – indices into eigvec for identifying states in eig_infected
  • all_drop_eigen_idx – indices into state and lin_state for identifying states not in eigvec
  • all_drop_infected_idx – indices into state and lin_state for identifying states not in eig_infected
  • eigen_drop_infected_idx – indices into eigvec for identifying states not in eig_infected
  • susceptible_idx – indices into state and lin_state for identifying states associated with the initial susceptible population
• initial_population – indices into the parameter vector for identifying the total number of individuals in the population and the initial total number of infected individuals in the population
  • total_idx
  • infected_idx

The justification for the initialization_mapping is given here. One potential source of confusion in the above process is the need to keep track of the different subsets of states. To help clarify this point, there are conceptually three types of state vectors containing nested subsets of states. 1. all_states – vectors containing all states 2. eigen_states – vectors containing all states of the linearized system that are included in the dominant eigenvector 3. infected_states – vectors containing all infected states, which are eigen_states that will be initialized using corresponding elements in the eigenvector

To make life a little easier for development on the C++ side we provide index vectors for converting between the three nested sets of states. 1. all_to_eigen_idx 2. all_to_infected_idx 3. eigen_to_infected_idx

And for completeness, in case it is useful, we also provide index vectors for dropping states from each kind of state vector. 1. all_drop_eigen_idx 2. all_drop_infected_idx 3. eigen_drop_infected_idx

struc objects

Sums of state variables and parameters

Within each simulation iteration, a user-defined number of sums of state variables and parameters are computed and stored immediately before the rate matrix is updated.

struc objects

As an example, the force of infection under the two-dose vaccination model can be given by the following struc objects.

Start by getting a set of vaccination parameters and state variables.

Create a symbolic vector of I-state variables, with 4 base states by 5 vaccination categories = 20 variables.

base_params <- read_params("PHAC.csv")
vax_params <- expand_params_vax(
  params = base_params,
  model_type = "twodose"
)
Istate = (McMasterPandemic:::expand_names(
  c('Ia', 'Ip', 'Im', 'Is'),   # = Icats
  attr(vax_params, 'vax_cat')) # = vax_cats
  %>% struc
)
as.matrix(Istate)

##       [,1]              
##  [1,] "(Ia_unvax)"      
##  [2,] "(Ip_unvax)"      
##  [3,] "(Im_unvax)"      
##  [4,] "(Is_unvax)"      
##  [5,] "(Ia_vaxdose1)"   
##  [6,] "(Ip_vaxdose1)"   
##  [7,] "(Im_vaxdose1)"   
##  [8,] "(Is_vaxdose1)"   
##  [9,] "(Ia_vaxprotect1)"
## [10,] "(Ip_vaxprotect1)"
## [11,] "(Im_vaxprotect1)"
## [12,] "(Is_vaxprotect1)"
## [13,] "(Ia_vaxdose2)"   
## [14,] "(Ip_vaxdose2)"   
## [15,] "(Im_vaxdose2)"   
## [16,] "(Is_vaxdose2)"   
## [17,] "(Ia_vaxprotect2)"
## [18,] "(Ip_vaxprotect2)"
## [19,] "(Im_vaxprotect2)"
## [20,] "(Is_vaxprotect2)"

Create a vector of base transmission rates associated with each of the 4 base I-states.

baseline_trans_rates =
  struc(
  'Ca',
  'Cp',
  '(1 - iso_m) * (Cm)',
  '(1 - iso_s) * (Cs)') *
  struc('(beta0) * (1/N)')
as.matrix(baseline_trans_rates)

##      [,1]                                  
## [1,] "(Ca) * (beta0) * (1/N)"              
## [2,] "(Cp) * (beta0) * (1/N)"              
## [3,] "(1 - iso_m) * (Cm) * (beta0) * (1/N)"
## [4,] "(1 - iso_s) * (Cs) * (beta0) * (1/N)"

Create a 5-by-5 block matrix of vaccine efficacies

vax_trans_red = struc_block(struc(
  '1',
  '1',
  '(1 - vax_efficacy_dose1)',
  '(1 - vax_efficacy_dose1)',
  '(1 - vax_efficacy_dose2)'),
  row_times = 1, col_times = 5)
as.matrix(vax_trans_red)

##      [,1]                       [,2]                      
## [1,] "(1)"                      "(1)"                     
## [2,] "(1)"                      "(1)"                     
## [3,] "(1 - vax_efficacy_dose1)" "(1 - vax_efficacy_dose1)"
## [4,] "(1 - vax_efficacy_dose1)" "(1 - vax_efficacy_dose1)"
## [5,] "(1 - vax_efficacy_dose2)" "(1 - vax_efficacy_dose2)"
##      [,3]                       [,4]                      
## [1,] "(1)"                      "(1)"                     
## [2,] "(1)"                      "(1)"                     
## [3,] "(1 - vax_efficacy_dose1)" "(1 - vax_efficacy_dose1)"
## [4,] "(1 - vax_efficacy_dose1)" "(1 - vax_efficacy_dose1)"
## [5,] "(1 - vax_efficacy_dose2)" "(1 - vax_efficacy_dose2)"
##      [,5]                      
## [1,] "(1)"                     
## [2,] "(1)"                     
## [3,] "(1 - vax_efficacy_dose1)"
## [4,] "(1 - vax_efficacy_dose1)"
## [5,] "(1 - vax_efficacy_dose2)"

And finally the force of infection vector, which looks like we really need to introduce the capability for common factors.

foi = kronecker(vax_trans_red, t(baseline_trans_rates)) %*% Istate
foi

## struc object with 5 rows and 1 column:
## 
## Row 1 
## (Ca) * (beta0) * (1/N) * (Ia_unvax) +
##  (Cp) * (beta0) * (1/N) * (Ip_unvax) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_unvax) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_unvax) +
##  (Ca) * (beta0) * (1/N) * (Ia_vaxdose1) +
##  (Cp) * (beta0) * (1/N) * (Ip_vaxdose1) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose1) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose1) +
##  (Ca) * (beta0) * (1/N) * (Ia_vaxprotect1) +
##  (Cp) * (beta0) * (1/N) * (Ip_vaxprotect1) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect1) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect1) +
##  (Ca) * (beta0) * (1/N) * (Ia_vaxdose2) +
##  (Cp) * (beta0) * (1/N) * (Ip_vaxdose2) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose2) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose2) +
##  (Ca) * (beta0) * (1/N) * (Ia_vaxprotect2) +
##  (Cp) * (beta0) * (1/N) * (Ip_vaxprotect2) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect2) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect2)
## Row 2 
## (Ca) * (beta0) * (1/N) * (Ia_unvax) +
##  (Cp) * (beta0) * (1/N) * (Ip_unvax) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_unvax) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_unvax) +
##  (Ca) * (beta0) * (1/N) * (Ia_vaxdose1) +
##  (Cp) * (beta0) * (1/N) * (Ip_vaxdose1) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose1) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose1) +
##  (Ca) * (beta0) * (1/N) * (Ia_vaxprotect1) +
##  (Cp) * (beta0) * (1/N) * (Ip_vaxprotect1) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect1) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect1) +
##  (Ca) * (beta0) * (1/N) * (Ia_vaxdose2) +
##  (Cp) * (beta0) * (1/N) * (Ip_vaxdose2) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose2) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose2) +
##  (Ca) * (beta0) * (1/N) * (Ia_vaxprotect2) +
##  (Cp) * (beta0) * (1/N) * (Ip_vaxprotect2) +
##  (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect2) +
##  (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect2)
## Row 3 
## (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_unvax) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_unvax) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_unvax) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_unvax) +
##  (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_vaxdose1) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_vaxdose1) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose1) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose1) +
##  (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_vaxprotect1) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_vaxprotect1) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect1) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect1) +
##  (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_vaxdose2) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_vaxdose2) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose2) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose2) +
##  (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_vaxprotect2) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_vaxprotect2) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect2) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect2)
## Row 4 
## (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_unvax) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_unvax) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_unvax) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_unvax) +
##  (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_vaxdose1) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_vaxdose1) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose1) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose1) +
##  (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_vaxprotect1) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_vaxprotect1) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect1) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect1) +
##  (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_vaxdose2) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_vaxdose2) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose2) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose2) +
##  (1 - vax_efficacy_dose1) * (Ca) * (beta0) * (1/N) * (Ia_vaxprotect2) +
##  (1 - vax_efficacy_dose1) * (Cp) * (beta0) * (1/N) * (Ip_vaxprotect2) +
##  (1 - vax_efficacy_dose1) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect2) +
##  (1 - vax_efficacy_dose1) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect2)
## Row 5 
## (1 - vax_efficacy_dose2) * (Ca) * (beta0) * (1/N) * (Ia_unvax) +
##  (1 - vax_efficacy_dose2) * (Cp) * (beta0) * (1/N) * (Ip_unvax) +
##  (1 - vax_efficacy_dose2) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_unvax) +
##  (1 - vax_efficacy_dose2) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_unvax) +
##  (1 - vax_efficacy_dose2) * (Ca) * (beta0) * (1/N) * (Ia_vaxdose1) +
##  (1 - vax_efficacy_dose2) * (Cp) * (beta0) * (1/N) * (Ip_vaxdose1) +
##  (1 - vax_efficacy_dose2) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose1) +
##  (1 - vax_efficacy_dose2) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose1) +
##  (1 - vax_efficacy_dose2) * (Ca) * (beta0) * (1/N) * (Ia_vaxprotect1) +
##  (1 - vax_efficacy_dose2) * (Cp) * (beta0) * (1/N) * (Ip_vaxprotect1) +
##  (1 - vax_efficacy_dose2) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect1) +
##  (1 - vax_efficacy_dose2) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect1) +
##  (1 - vax_efficacy_dose2) * (Ca) * (beta0) * (1/N) * (Ia_vaxdose2) +
##  (1 - vax_efficacy_dose2) * (Cp) * (beta0) * (1/N) * (Ip_vaxdose2) +
##  (1 - vax_efficacy_dose2) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxdose2) +
##  (1 - vax_efficacy_dose2) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxdose2) +
##  (1 - vax_efficacy_dose2) * (Ca) * (beta0) * (1/N) * (Ia_vaxprotect2) +
##  (1 - vax_efficacy_dose2) * (Cp) * (beta0) * (1/N) * (Ip_vaxprotect2) +
##  (1 - vax_efficacy_dose2) * (1 - iso_m) * (Cm) * (beta0) * (1/N) * (Im_vaxprotect2) +
##  (1 - vax_efficacy_dose2) * (1 - iso_s) * (Cs) * (beta0) * (1/N) * (Is_vaxprotect2)

Sums of state variables and parameters

Start with a basic example model.

params <- read_params("ICU1.csv")
state <- make_state(params = params)
M <- McMasterPandemic::make_ratemat(state, params, sparse = TRUE)
test_model <- flexmodel(
    params, state,
    start_date = "2021-09-09", end_date = "2021-10-09"
)

Then we define sums of state variables and parameter vectors. Here we add the sum of the I states and of the ICU states using regular expressions and just a list of state names.

test_model = (test_model
  %>% add_state_param_sum(sum = "Isum", summands = "^I[a-z]")
  %>% add_state_param_sum(sum = "ICUsum", summands = c("ICUs", "ICUd"))
)

Then we can just proceed as we normally would, but with the option to refer to these sums (this is not meant to reflect real epidemiology, but to illustrate the interface).

options(MP_flex_spec_version = "0.0.2")
test_model = (test_model
  %>% add_rate("E", "Ia", ~ (alpha) * (sigma) * (1 / Isum))
  %>% add_rate("E", "Ip", ~ (1 - alpha) * (sigma))
  %>% add_rate("Ia", "R", ~ (gamma_a))
  %>% add_rate("Ip", "Im", ~ (mu) * (gamma_p))
  %>% add_rate("Ip", "Is", ~ (1 - mu) * (gamma_p))
  %>% add_rate("Im", "R", ~ (gamma_m))
  %>% add_rate("Is", "H", ~ (1 - nonhosp_mort) * (phi1) * (gamma_s))
  %>% add_rate("Is", "ICUs", ~ (1 - nonhosp_mort) * (1 - phi1) * (1 - phi2) * (gamma_s))
  %>% add_rate("Is", "ICUd", ~ (1 - nonhosp_mort) * (1 - phi1) * (phi2) * (gamma_s))
  %>% add_rate("Is", "D", ~ (nonhosp_mort) * (gamma_s))
  %>% add_rate("ICUs", "H2", ~ (psi1))
  %>% add_rate("ICUd", "D", ~ (psi2))
  %>% add_rate("H2", "R", ~ (psi3))
  %>% add_rate("H", "R", ~ (rho) / (1 / ICUsum))
  %>% add_rate("Is", "X", ~ (1 - nonhosp_mort) * (phi1) * (gamma_s))
  %>% add_rate("S", "E", ~
                 (Ia) * (beta0) * (1 / N) * (Ca) +
                 (Ip) * (beta0) * (1 / N) * (Cp) +
                 (Im) * (beta0) * (1 / N) * (Cm) * (1 - iso_m) +
                 (Is) * (beta0) * (1 / N) * (Cs) * (1 - iso_s))
  %>% add_parallel_accumulators(c("X", "N", "P", "V"))
  %>% update_tmb_indices()
)

And here are the indices that are relevant to the tracking of sums of state variables and parameters.

test_model$tmb_indices$sum_indices

## $sumidx
## integer(0)
## 
## $sumcount
## integer(0)
## 
## $summandidx
## integer(0)

We describe these below.

struc objects

Sums of state variables and parameters

We introduce three new index/count vectors.

• sumidx – indices, locating the variable to store each sum, into the state_param vector defined in version 0.0.1
• sumcount – vector containing the number of summands required for each sum
• summandidx – indices, locating the summands, into the state_param vector