Package 'macpan2'

Description

Convert a function that represents an elementwise binary operator into one that is consistent with the C++ engine. This function is intended to clarify how macpan2 treats binary operators, which is a little different from base R. The difference is clarified in vignette("elementwise_binary_operators"), and BinaryOperator is primarily used as a resource for that vignette.

Usage

BinaryOperator(operator)

mp_binary_operator(operator)

Arguments

Value

A binary operator consistent with the C++ engine.

Functions

• BinaryOperator(): Synonym of mp_binary_operator.

Examples

set.seed(1L)
A = matrix(abs(rnorm(6)), 3, 2)  # 3 by 2 matrix
x = matrix(abs(rnorm(3)))        # 3 by 1 matrix
y = t(abs(rnorm(2)))             # 1 by 2 matrix
times = BinaryOperator(`*`)
pow = BinaryOperator(`^`)
identical(times(A, x), times(x, A))  ## TRUE
identical(pow(A, y), pow(y, A))  ## FALSE

Description

Comparison Functions

Usage

all_equal(x, y)

all_consistent(x, y)

not_all_equal(x, y)

all_not_equal(x, y)

Arguments

Functions

• all_equal(): Is it true that all corresponding elements of x and y are equal, have the same shape, and have no missing values?

• all_consistent(): Is it true that all corresponding elements of x and y are either equal or at least one is a blank string, have the same shape, and have no missing values?

• not_all_equal(): Complement of all_equal.

• all_not_equal(): Do not know yet. Currently unused; should we remove?

Description

Distributions which can be used to specify prior or likelihood components in model calibration.

• Uniform Distribution (Improper), only appropriate for prior components - mp_unif

• Normal Distribution - mp_norm

• Log-Normal Distribution - mp_log_norm

• Logit-Normal Distribution - mp_logitnorm

• Poisson Distribution - mp_pois

• Negative Binomial Distribution - mp_nbinom

Usage

mp_unif(trans_distr_param = list())

mp_uniform(trans_distr_param = list())

mp_norm(
  location = mp_distr_param_null("location"),
  sd,
  trans_distr_param = list(location = mp_identity, sd = mp_log)
)

mp_normal(
  location = mp_distr_param_null("location"),
  sd,
  trans_distr_param = list(location = mp_identity, sd = mp_log)
)

mp_lnorm(
  location = mp_distr_param_null("location"),
  sd,
  trans_distr_param = list(location = mp_identity, sd = mp_identity)
)

mp_log_normal(
  location = mp_distr_param_null("location"),
  sd,
  trans_distr_param = list(location = mp_identity, sd = mp_identity)
)

mp_logitnorm(
  location = mp_distr_param_null("location"),
  sd,
  trans_distr_param = list(location = mp_identity, sd = mp_identity)
)

mp_logit_normal(
  location = mp_distr_param_null("location"),
  sd,
  trans_distr_param = list(location = mp_identity, sd = mp_identity)
)

mp_pois(
  location = mp_distr_param_null("location"),
  trans_distr_param = list(location = mp_identity)
)

mp_poisson(
  location = mp_distr_param_null("location"),
  trans_distr_param = list(location = mp_identity)
)

mp_nbinom(
  location = mp_distr_param_null("location"),
  disp,
  trans_distr_param = list(location = mp_identity, disp = mp_log)
)

mp_neg_bin(
  location = mp_distr_param_null("location"),
  disp,
  trans_distr_param = list(location = mp_identity, disp = mp_log)
)

Arguments

Details

All distributional parameter arguments can be specified either as a numeric value, a character string giving the parameter name, or a distributional parameter object (See fit_distr_params).

Description

Convert a compartmental model specification into a graph (using the graph package) of the box diagram for the model.

Usage

mp_dot_layout(
  spec,
  include_inout = FALSE,
  action = c("render", "layout", "define"),
  define_args = list(edgemode = "directed"),
  layout_args = list(attrs = list(graph = list(rankdir = "LR"), node = list(shape =
    "rectangle"))),
  render_args = list(edgemode = "directed")
)

dot_layout(spec, include_inout = FALSE)

Arguments

Details

In order to create these graph objects you need to install the graph package, and to plot them on the current graphics device you need to install the Rgraphviz package.

Value

A graphAM object.

Functions

• dot_layout(): Deprecated in favour of mp_dot_layout, which both plots and returns the graphAM object.

Examples

("macpan_base"
  |> mp_official_library()
  |> mp_dot_layout()
)

Description

Empty matrices are useful when defining matrices that do not need to be initialized because they will get computed before they are required by other expressions. They can also provide a useful placeholder for matrices that should only have a value after a certain phase in the simulation.

Usage

empty_matrix

Format

A numeric matrix with zero rows and zero columns.

Examples

spec = mp_tmb_model_spec(during = list(x ~ time_step(0)))
identical(spec$empty_matrices()$x, empty_matrix) ## TRUE

Description

Output of mp_trajectory if nothing is simulated.

Usage

empty_trajectory

Format

A data frame with zero rows and the following columns: matrix, time, row, col, value.

Description

Evaluate an expression in the TMB-based C++ simulation and objective function engine. This function is useful for trying out the engine_functions that can be used to define macpan2 models.

Usage

engine_eval(
  .expr,
  ...,
  .matrix_to_return,
  .tmb_cpp = getOption("macpan2_dll"),
  .structure_labels = NullLabels()
)

Arguments

Value

Matrix being produced on the right-hand-side or matrix given in .matrix_to_return if it is provided.

Examples

engine_eval(~ exp(y), y = pi) # ~ 23.14069

# It is not currently possible to assign values to a subset of
# a matrix in a natural way (e.g. you cannot do things like x[1] = exp(y)),
# but you can achieve this functionality using the assign function.
engine_eval(~ assign(x, 1, 0, exp(y))
  , x = rep(0, 2)
  , y = pi
  , .matrix_to_return = "x"
)

Description

Functions currently supported by the C++ TMB engine for constructing expressions for defining model simulations.

Details

The quickest way to experiment with these functions is to use the engine_eval function, as in the following example that calculates a force of infection.

engine_eval(~ beta * I / N
  , beta = 0.25
  , I = 1e3
  , N = 1e7
)

To produce a dynamical simulation that iteratively evaluates expression involving these functions, use simple_sims.

simple_sims(
  iteration_exprs = list(x ~ x - 0.9 * x),
  time_steps = 5,
  mats = list(x = 1)
)

Here, x - 0.9 * x is assigned to x at each of five iterations of a simulation loop.

If these expressions involve matrices with more than one element, You can control which elements in the evaluation of the right hand side go to which elements on the left hand side. This technique involves using either square brackets or the c function on the left hand side. For more information on assignment, please see the section on Assignment below.

Elementwise Binary Operators

Elementwise binary operators take two matrix-valued arguments and apply a binary operator (e.g. +, *) to each set of corresponding elements, and return the corresponding matrix-valued output containing the resulting elements. What does 'corresponding' mean? If the two matrix-valued arguments have the same shape (same number of rows and columns), then two elements correspond if they occur in the same row and column position in the two matrices. If the two matrices are not of the same shape but there is one row and/or one column in either matrix, then the singleton rows and columns are recycled sufficiently many times so that they match the shape of the other matrix. If after recycling singleton rows and columns the matrices are still of different shape, then an error is thrown and the matrices are said to be incompatible.

Functions

• x + y

• x - y

• x * y

• x / y

• x ^ y

Arguments

• x : Any matrix with dimensions compatible with y.

• y : Any matrix with dimensions compatible with x.

Return

• A matrix with the binary operator applied elementwise after any necessary recycling of rows and/or columns.

Examples

engine_eval(~ 1 + 2)
engine_eval(~ y * z, y = 1:3, z = matrix(1:6, 3, 2))
engine_eval(~ 1 / (1 - y), y = 1/4)

Elementwise Math

Functions

• log(x) : Natural logarithm.

• exp(x) : Exponential function.

• cos(x) : Cosine function.

• sin(x) : Sine function.

• sqrt(x) : Squareroot function.

• invlogit(x) : Inverse logit function, 1/(1 + exp(-x)).

• logit(x) : Logit function, log(x/(1-x)).

Arguments

• x : Any numeric matrix.

Return

• A matrix with the same dimensions as x, containing the results of applying the function to each element of x.

Examples

engine_eval(~ log(y), y = c(2, 0.5))

Limiting Values

Functions

• proportions(x, limit, eps)

• limit(x, limit, eps)

• clamp(x, eps, limit)

Arguments

• x : Any matrix

• limit : numeric value to return elementwise from proportions if sum(x) < eps

• eps : numeric tolerance for sum(x)

Return

• matrix of x / sum(x) or rep(limit, length(x)) if sum(x) < eps.

Examples

engine_eval(~ proportions(y, 0.5, 1e-8), y = c(2, 0.5))

Details

The return value depends on the function.

• proportions : matrix of x / sum(x) or rep(limit, length(x)) if sum(x) < eps.

• limit :

• clamp :

Details

The divide_safe() function conducts elementwise division of the first two arguments, and then allows two other arguments:

• limit : Value to return if the sqrt(denominator^2) is

Details

The clamp function smoothly clamps the elements of a matrix so that they do not get closer to 0 than a tolerance, eps, with a default of 1e-12. This clamp function is the following modification of the squareplus function.

$f(x) = \epsilon_- + \frac{(x - \epsilon_-) + \sqrt{(x - \epsilon_-)^2 + (2\epsilon_0 - \epsilon_-)^2 - \epsilon_-^2}}{2}$

Where the two parameters are defined as follows.

$\epsilon_0 = f(0)$

$\epsilon_- = \lim_{x \to -\infty}f(x)$

This function is differentiable everywhere, monotonically increasing, and $f(x) \approx x$ if $x$ is positive and not too close to zero. By modifying the parameters, you can control the distance between $f(x)$ and the horizontal axis at two 'places' – $0$ and $-\infty$. See issue #93. for more information.

For clamp the arguments specifically mean.

• x : A matrix with elements that should remain positive.

• eps : A small positive number, $\epsilon_0 = f(0)$, giving the value of the function when the input is zero. The default value is 1e-11

• limit : A small positive number,

  $\epsilon_- = \lim_{x \to -\infty}f(x)$

  , giving the value of the function as the input goes to negative infinity. The default is limit = 1e-12. This limit should be chosen to be less than eps to ensure that clamp is twice differentiable.

Integer Sequences

Functions

• from:to : Inclusive and ordered sequence of integers between two bounds.

• seq(from, length, by) : Ordered sequence of integers with equal spacing between adjacent values.

Arguments

• from : Scalar integer giving the first integer in the sequence.

• to : Scalar integer giving the last integer in the sequence.

• length : Number of integers in the sequence.

• by : Scalar giving the difference between adjacent values in the sequence.

Return

• Column vector with a sequence of integers.

Details

The colon operator works much like the base R version :. It takes two scalar-valued integers and returns a column vector with all integers between the two inputs.

The seq function is a little different from the base R default, seq, in that it allows the user precise control over the length of the output through the length argument. The base R function gives the user this option, but not as the default.

Examples

engine_eval(~ 1:10)
engine_eval(~ seq(1, 10, 2))

Replicate Elements

Functions

• rep(x, times) : Replicate a column vector a number of times, by repeatedly stacking it on top of itself.

• recycle(x, rows, cols) : Recycle rows and columns of x to produce a matrix with rows rows and cols columns.

Arguments

• x : A scalar-valued variable to repeat.

• times : A scalar-valued integer variable giving the number of times to repeat x.

• rows : Number of rows in the output of recycle.

• cols : Number of columns in the output of recycle.

Return

• Column vector with times copies of x stacked on top of each other.

Examples

engine_eval(~ rep(1, 10))
engine_eval(~ recycle(  1:3,  3, 4))
engine_eval(~ recycle(t(1:4), 3, 4))

Matrix Multiplication

Functions

• x %*% y : Standard matrix multiplication.

• x %x% y : Kronecker product

• sparse_mat_mult(x, i, j, y, z) : Matrix multiplication when the left matrix is represented as a column vector, x, of non-zero elements and integer vectors of row, i, and column, j, indices. The right matrix and the resulting matrix are not represented as sparse matrices.

Arguments

• x : A matrix.

• y : A matrix.

• i : Integer vector the same length as x giving zero-based row indices for sparse matrix representation.

• j : Integer vector the same length as x giving zero-based column indices for sparse matrix representation.

• z : A matrix with dimensions equal to the result of the sparse matrix multiplication (see details).

Return

• The matrix product of x and y.

Examples

engine_eval(~ (1:10) %*% t(1:10))
engine_eval(~ (1:10) %x% t(1:10))

Details

For standard matrix multiplication, x %*% y, the number of columns of x equals the number of rows of y.

Think about sparse_mat_mult(x, i, j, y, z) as similar to z ~ x %*% y, where x is represented differently. In particular, the argument x is a column vector containing the non-zero elements of the left matrix, i contains the zero-based row indices associated with each element in x, and j contains the zero-based column indices associated with each element in x.

Parenthesis

The order of operations can be enforced in the usual way with round parentheses, (.

Reshaping and Combining Matrices

Functions

• c(...) : Stack columns of arguments into a single column vector.

• cbind(...) : Create a matrix containing all of the columns of a group of matrices with the same number of rows.

• rbind(...) : Create a matrix containing all of the rows of a group of matrices with the same number of columns.

• matrix(x, rows, cols) : Reshape a matrix to have rows rows and cols columns. The input x must have rows * cols elements.

• t(x) : Standard matrix transpose.

Arguments

• ... : Any number of dimensionally consistent matrices. The definition of dimensionally consistent depends on the function.

• x : Can be any matrix for t, but for matrix it must have rows * cols elements.

• rows : Scalar integer giving the number of rows in the output.

• cols : Scalar integer giving the number of columns in the output.

Return

• A combined or reshaped matrix.

Details

Any number of column vectors can be combined into a bigger column vector.

Column and row vectors of the same length can be combined using the cbind and rbind functions respectively

The matrix function can be used to redefine the numbers of rows and columns to use for arranging the values of a matrix. It works similarly to the base R matrix function in that it takes the same arguments. On the other hand, this function differs substantially from the base R version in that it must be filled by column and there is no byrow option.

Matrices can be transposed with the usual function, t.

Examples

engine_eval(~ c(a, b, c), a = 1, b = 10:13, c = matrix(20:25, 3, 2))
engine_eval(~ cbind(a, 10 + a), a = 0:3)
engine_eval(~ rbind(a, 10 + a), a = t(0:3))
engine_eval(~ matrix(1:12, 4, 3))
engine_eval(~ t(1:3))

Matrix Diagonals

Functions

• to_diag(x) : Create a diagonal matrix by setting the diagonal to a column vector, x.

• from_diag(x) : Extract the diagonal from a matrix, x, and return the diagonal as a column vector.

Arguments

• x : Any matrix (for from_diag) or a column vector (for to_diag). It is common to assume that x is square for from_diag but this is not required.

Return

• to_diag(x) : Diagonal matrix with x on the diagonal.

• from_diag(x) : Column vector containing the diagonal of x. A value is considered to be on the diagonal if it has a row index equal to the column index.

Details

The to_diag function can be used to produce a diagonal matrix by setting a column vector equal to the desired diagonal. The from_diag does (almost) the opposite, which is to get a column vector containing the diagonal of an existing matrix.

Examples

engine_eval(~from_diag(matrix(1:9, 3, 3)))
engine_eval(~to_diag(from_diag(matrix(1:9, 3, 3))))
engine_eval(~from_diag(to_diag(from_diag(matrix(1:9, 3, 3)))))

Summarizing Matrix Values

Functions

• sum(...) : Sum all of the elements of all of the matrices passed to ....

• col_sums(x) : Row vector containing the sums of each column.

• row_sums(x) : Column vector containing the sums of each row.

• group_sums(x, f, n) : Column vector containing the sums of groups of elements in x. The groups are determined by the integers in f and the order of the sums in the output is determined by these integers.

• mean(x) : Arthmetic average of all elements in matrix x.

• sd(x) : Sample standard deviation of all elements in matrix x.

Arguments

• ... : Any number of matrices of any shape.

• x : A matrix of any dimensions, except for group_sums that expects x to be a column vector.

• f : A column vector the same length as x containing integers between 0 and m-1, given m unique groups. Elements of f refer to the indices of x that will be grouped and summed.

• n : A column vector of length m. If f does not contain group k in ⁠[0, m-1]⁠, group_sums skips this group and the output at index k+1 is n[k+1].

Return

• A matrix containing sums of subsets of the inputs.

Details

The row_sums and col_sums are similar to the base R rowSums and colSums functions, but with slightly different behaviour. In particular, the row_sums function returns a column vector and the col_sums function returns a row vector. If a specific shape is required then the transpose t function must be explicitly used.

Examples

x = 1
y = 1:3
A = matrix(1:12, 4, 3)
engine_eval(~ sum(y), y = y)
engine_eval(~ sum(x, y, A), x = x, y = y, A = A)
engine_eval(~ col_sums(A), A = A)
engine_eval(~ row_sums(A), A = A)
engine_eval(~ group_sums(x, f, n), x = 1:10, f = rep(0:3, 1:4), n = c(1:4))

Sweeping Matrix Elements

Functions

• cumsum(x) : Return a matrix with columns containing the cumulative sum of the columns in x.

Arguments

• x : A matrix.

Return

A matrix the same size as x but with columns containing the cumulative sum of the columns in x.

Extracting Matrix Elements

Functions

• x[i,j] : Return a matrix containing a subset of the rows and columns of x.

• block(x,i,j,n,m) : Return a matrix containing a contiguous subset of rows and columns of x https://eigen.tuxfamily.org/dox/group__TutorialBlockOperations.html.

• last(x) : The last element of a matrix (i.e., the lower-right element).

Arguments

• x : Any matrix.

• i : An integer column vector (for [) or integer scalar (for block) containing the indices of the rows to extract (for [) or the index of the first row to extract (for block).

• j : An integer column vector (for [) or integer scalar (for block) containing the indices of the columns to extract (for [) or the index of the first column to extract (for block). If j is missing in a call to [, it is assumed to be j = 0 although we might change this default to be the vector of all column indices.

• n : Number of rows in the block to return.

• m : Number of columns in the block to return.

Return

• A matrix containing a subset of the rows and columns in x.

Details

Note that zero-based indexing is used so the first row/column gets index, 0, etc. The block function is expected to be more efficient than [ when the elements to be extracted are contiguous.

Examples

engine_eval(~ A[c(3, 1, 2), 2], A = matrix(1:12, 4, 3))
engine_eval(~ block(x,i,j,n,m), x = matrix(1:12, 4, 3), i=1, j=1, n=2, m=2)
engine_eval(~ last(A), A = matrix(1:12, 4, 3))

Accessing Past Values in the Simulation History

For matrices with their simulation history saved, it is possible to bind the rows or columns of past versions of such matrices into a single matrix.

Functions

• rbind_lag(x, lag, t_min) : Bind the rows of versions of x that were recorded at the end of all simulation iterations corresponding to time lags given by integers in lag.

• rbind_time(x, t, t_min) : Bind the rows of versions of x that were recorded at the end of all simulation iterations corresponding to integers in t.

• cbind_lag(x, lag, t_min) : Bind the columns of versions of x that were recorded at the end of all simulation iterations corresponding to time lags given by integers in lag. (TODO – cbind_lag is not developed yet)

• cbind_time(x, t, t_min) : Bind the columns of versions of x that were recorded at the end of all simulation iterations corresponding to integers in t. (TODO – cbind_lag is not developed yet)

Arguments

• x : Any matrix with saved history such that the number of columns (for ⁠rbind_*⁠) or rows (for ⁠cbind_*⁠) does not change throughout the simulation.

• lag : Integer vector giving numbers of time steps before the current step to obtain past values of x.

• t : Integer vector giving time steps at which to obtain past values of x.

• t_min : Integer giving the minimum time step that is allowed to be accessed. All time-steps in t or implied by lag that are before t_min are ignored.

Return

• A matrix containing values of x from past times.

Time Indexing

Get or update the index of the current or lagged time step or the index of the current time group. A time group is a contiguous set of time steps defined by two change points.

Functions

• time_step(lag): Get the time-step associated with a particular lag from the current time-step. If the lagged time-step is less than zero, the function returns zero.

• time_group(index, change_points): Update the index associated with the current time group. The current group is defined by the minimum of all elements of change_points that are greater than the current time step. The time group index is the index associated with this element. Please see the examples below, they are easier to understand than this explanation.

• time_var(x, change_points): An improvement to time_group. Returns values in x at time steps in change_points, return value remains constant between change_points.

Arguments

• x: Column vector representing a time series. time_var will return the value of x corresponding to element in change_points that contains the current time.

• lag: Number of time-steps to look back for the time-step to return.

• change_points: Increasing column vector of time steps giving the lower bound of each time group.

Return

A 1-by-1 matrix with the time-step lag steps ago, or with zero if t+1 < lag

Examples

simple_sims(
  iteration_exprs = list(x ~ time_step(0)),
  time_steps = 10,
  mats = list(x = empty_matrix)
)
sims = simple_sims(
  iteration_exprs = list(
    j ~ time_group(j, change_points),
    time_varying_parameter ~ time_variation_schedule[j]
  ),
  mats = list(
    j = 0,
    change_points = c(0, 4, 7),
    time_variation_schedule = c(42, pi, sqrt(2)),
    time_varying_parameter = empty_matrix
  ),
  time_steps = 10,
)
set.seed(1L)
change_points = c(0,2,5)
x_val = rnorm(length(change_points))
simple_sims(
    iteration_exprs = list(x ~ time_var(x_val,change_points))
  , int_vecs = list(change_points = change_points)
  , mats = list(x = empty_matrix, x_val=x_val)
  , time_steps = 10
)

Convolution

You can take the convolution of each element in a matrix, x, over simulation time using a kernel, k.

Functions

• convolution(x, k)

Arguments

• x : The matrix containing elements to be convolved.

• k : A column vector giving the convolution kernel.

Return

A matrix the same size as x but with the convolutions, $y_{ij}$, of each element, $x_{ij}$. The value of $y_{ij}$ at time $t = 1, ..., T$ is given by the following.

$y_{ij}(t) = \sum_{\tau = 0}^{min(t,m)-1} x_{ij}(t-\tau) k_\tau$

Where:

• $x_{ij}(t)$ : value of $x_{ij}$ at time step $t$.

• $y_{ij}(t)$ : value of $y_{ij}$ at time step $t$.

• $t = 1, ..., T$ : the time step.

• $\tau = 0, ..., m - 1$ : index of the time lag for a kernel of length $m$.

• $k_\tau$ : value of the kernel associated with lag $\tau$.

Details

If any empty matrices are encountered when looking back in time, they are treated as matrices with all zeros. The convolution of a matrix of all positive values will be biased low for all time steps less than the length of the kernel (i.e., for all time steps such that $t < m$), and therefore one should only compare observed data with a convolution (e.g., when calibrating) for time steps less than $m$.

Examples

simple_sims(
  list(
    x ~ 3 * x * (1 - x),
    y ~ convolution(x, rep(1/10, 10))
  ),
  time_steps = 50,
  mats = list(x = 0.5, y = empty_matrix)
)

Probability Densities

All probability densities have the same first two arguments.

• observed

• simulated

The simulated argument gives a matrix of means for the observed values at which the densities are being evaluated. Additional arguments are other distributional parameters such as the standard deviation or dispersion parameter. All densities are given as log-densities, so if you would like the density itself you must pass the result through the exp function.

If the simulated matrix or the additional parameter matrices have either a single row or single column, these singleton rows and columns are repeated to match the number of rows and columns in the observed matrix. This feature allows one to do things like specify a single common mean for several values.

Functions

• dpois(observed, simulated) : Log of the Poisson density based on this dpois TMB function.

• dnbinom(observed, simulated, over_dispersion) : Log of the negative binomial density based on this dnbinom TMB function. To get the variance that this function requires we use this expression, simulated + simulated^2/over_dispersion, following p.165 in this book

• dnorm(observed, simulated, standard_deviation) : Log of the normal density based on this dnorm TMB function.

• dbinom(observed, size, probability) : Log of the binomial density based on the dbinom TMB function.

Arguments

• observed : Matrix of observed values at which the density is being evaluated.

• simulated : Matrix of distributional means, with singleton rows and columns recycled to match the numbers of rows and columns in observed.

• over_dispersion : Over-dispersion parameter given by (simulated/standard_deviation)^2 - simulated).

• standard_deviation : Standard deviation parameter.

• size : Number of Bernoulli trials.

• probability : Probability of a successful Bernoulli trial.

Pseudo-Random Number Generators

All random number generator functions have mean as the first argument. Subsequent arguments give additional distributional parameters. Singleton rows and columns in the matrices passed to the additional distributional parameters are recycled so that all arguments have the same number of rows and columns. All functions return a matrix the same shape as mean but with pseudo-random numbers deviating from each mean in the mean matrix.

Functions

• rpois(mean) : Pseudo-random Poisson distributed values.

• rnbinom(mean, over_dispersion) : Pseudo-random negative binomially distributed values.

• rnorm(mean, standard_deviation) : Pseudo-random normal values.

• rbinom(size, prob) : Pseudo-random binomial values.

• reulermultinom(size, rate, dt) : Pseudo-random Euler-multinomial values.

Arguments

• mean : Matrix of means about which to simulate pseudo-random variation.

• over_dispersion : Matrix of over-dispersion parameters given by (simulated/standard_deviation)^2 - simulated).

• standard_deviation : Matrix of standard deviation parameters.

• size : Matrix of numbers of trials.

• prob : Matrix of probabilities of success.

• rate : Matrix of rates, used to compute the probabilities in a multinomial distribution. The probability associated with the ith rate, r_i, is (1 - exp(-sum(r * dt))) * (r_i / r), where r is the sum of the rates. This is not a typical multinomial distribution in that if you sum these probabilities up you do not get 1 but rather (1 - exp(-sum(r * dt))). See details below for more on the Euler-multinomial distribution

• dt : Optional parameter specifying the length of the time step. See details below for more on the Euler-multinomial distribution.

Details

The Euler-multinomial distribution is used to model how many individuals move from one compartment to a set of other compartments in a single time step of length dt. The rate of moving to each of these compartments is characterized by the associated element in the rate matrix. The reason why the probabilities do not sum to 1, is that not all individuals have to change compartments in a time step.

Cumulative Distribution Functions

Lower-tail cumulative distribution functions.

Functions

• pgamma(q, shape, scale) : Cumulative distribution function of the Gamma distribution. This is a lite wrapper for the pgamma function in TMB.

• pnorm(q, mean, sd) : Cumulative distribution function of the normal distribution. This is a lite wrapper for the pnorm function in TMB.

Arguments

• q : Matrix of Quantiles.

• shape : Matrix of shape parameters of the Gamma distribution.

• scale : Matrix of scale parameters of the Gamma distribution.

• mean : Matrix of mean parameters of the normal distribution.

• sd : Matrix of standard deviation parameters of the normal distribution.

Rounding

Functions

round(x) : Round elements of a matrix to the nearest integer.

Arguments

• x : Matrix to be rounded.

Details

Be careful if you are using rounding in a model to be calibrated. Rounding will break differentiability of the objective function if x depends, either directly or indirectly, on parameters being calibrated. This will lead to incorrect gradients potentially being passed to an optimizer. To be safe, do not round in models being calibrated.

Debugging Instrumentation

Functions to use when you are trying to figure stuff out.

Functions

• print(x) : Print out the value of a matrix.

• check_finite(x) : Stop the simulations and return an error if x has any non-finite values.

Arguments

• x : Name of a matrix in the model.

Return

An empty_matrix.

Examples

simple_sims(
     list(dummy ~ print(x), x ~ x / 2)
   , time_steps = 10
   , mats = list(x = 2)
)
engine_eval(~ 1/0) ## returns Inf
engine_eval(~ check_finite(1/0)) ## returns nothing and throws an error

Assign (deprecated)

Assign values to a subset of the elements in a matrix.

Functions

• assign(x, i, j, v)

Arguments

• x : Matrix with elements that are to be updated by the values in v.

• i : Column vector of row indices pointing to the elements of x to be updated. These indices are paired with those in v. If the length of i does not equal that of v, then it must have a single index that gets paired with every element of v. Indices are zero-based, i=0 corresponds to the first row.

• j : Column vector of column indices pointing to the elements of x to be updated. These indices are paired with those in v. If the length of j does not equal that of v, then it must have a single index that gets paired with every element of v. Indices are zero-based, j=0 corresponds to the first column.

• v : Column vector of values to replace elements of x at locations given by i and j.

Return

The assign function is not called for its return value, which is an empty_matrix, but rather to modify x but replacing some of its components with those in v.

Examples

x = matrix(1:12, 3, 4)
engine_eval(~ x + 1, x = x)
engine_eval(~ x + 1, x = x, .matrix_to_return = "x")
engine_eval(~ assign(x, 2, 1, 100), x = x, .matrix_to_return = "x")
engine_eval(~ assign(x
  , c(2, 1, 0)
  , 0
  , c(100, 1000, 10000)
), x = x, .matrix_to_return = "x")

Unpack (deprecated)

Unpack elements of a matrix into smaller matrices.

Functions

• unpack(x, ...)

Arguments

• x : Matrix with elements to be distributed to the matrices passed through ....

• ... : Matrices with elements to be replaced by the values of elements in x in column-major order. These matrices must be named matrices and not computed on the fly using expressions. Note that even subsetting (e.g. unpack(x, y[0], y[3])) counts as an expression. This use-case would require the assign function assign(y, c(0, 3), 0, x).

Return

The unpack function is not called for its return value, which is an empty_matrix, but rather to modify the matrices in ... by replacing at least some of its components with those in x.

Examples

Here we fill a matrix with integers from 1 to 12 and then unpack them one-at-a-time into two column vectors, x and y. By returning y we see the integers after the first three were used up by x.

engine_eval(~unpack(matrix(1:12, 3, 4), x, y)
  , x = rep(0, 3)
  , y = rep(1, 5)
  , .matrix_to_return = "y"
)

Assignment

The left-hand-side of formulas sent to the simulation engine determine assignment works.

Functions

• y ~ x : Assign x to y.

• y[i] ~ x : Assign the first column of x to those rows in the first column of y that are indexed by i.

• y[i, j] ~ x : Assign each element, x[k, l], in x, to element, y[i[k], j[l]], in y.

• c(...) ~ x : Assign the elements of the columns of x (stacked on top of each other) to the matrices in ... in the order in which they appear. If the number of columns is x equals the number of matrices in ..., and if these matrices are vectors (i.e., have only a single column or a single row), then the columns of x become assigned to the vectors in ....

Arguments

• x : Matrix containing the result of the expression on the right-hand-side.

• y : Matrix with elements that will be assigned the elements of x.

• i : Integer vector giving zero-based row indexes describing the rows in x that get the

Description

Find all paths through a compartmental model.

Usage

find_all_paths(edges_df, start_node_guesses = character(0L))

Arguments

Value

List of character vectors of state variable names, each vector describing a path through the model.

Examples

spec = mp_tmb_library("starter_models", "macpan_base", package = "macpan2")
spec |> mp_flow_frame() |> find_all_paths()

Description

Distributional parameters can be added to the list of parameters that are fit during calibration. By default, distributional parameters in priors and likelihoods are not fit. Use mp_nofit to exclude distributional parameters from being fit and mp_fit to fit them.

Usage

mp_fit(x, trans = DistrParamTransDefault())

mp_nofit(x, trans = DistrParamTransDefault())

Arguments

Value

A distributional parameter object.

Examples

# First we call the SIR model spec, and generate some data for calibration.
spec = mp_tmb_library("starter_models", "sir", package = "macpan2")
data = mp_simulator(spec, 50, "infection") |> mp_trajectory()

# Suppose we want to specify a Normal prior on the transmission parameter 
# beta, and we are interested in estimating the prior standard deviation.
# Here we use `mp_fit` to estimate the standard deviation, `sd`, and we 
# provide a numeric starting value for `sd` in the optimization. 
cal = mp_tmb_calibrator(
    spec
  , data
  , traj = "infection"
  , par = list(beta = mp_normal(location = 0.35, sd = mp_fit(0.1)))
  , default = list(beta = 0.25)
)

# When viewing the calibration objective function we can see the additional
# prior density term added for beta. The standard deviation parameter has
# been automatically named 'distr_params_log_sd_beta'.
cal$simulator$tmb_model$obj_fn$obj_fn_expr

# Next we optimize and view the fitted parameters. We can see the 
# distributional parameter in the coefficient table with a default value 
# equal to the numeric value we provided to `mp_fit` above.
mp_optimize(cal)
if (suppressPackageStartupMessages(require(broom.mixed))) {
  print(mp_tmb_coef(cal))
}

# If instead we want control over the name of the new fitted distributional
# parameter, we can add a new variable to our model specification with the 
# default value set to the desired optimization starting value.
updated_spec = spec |> mp_tmb_insert(default = list(sd_var = 0.1))

# In the calibrator, we use the name of this newly added variable, "sd_var",
# as input to `mp_fit`.
cal = mp_tmb_calibrator(
    updated_spec
  , data
  , traj = "infection"
  , par = list(beta = mp_normal(location = 0.35, sd = mp_fit("sd_var")))
  , default = list(beta = 0.25)
)

# We can see this distributional parameter get propogated to the objective 
# function and the fitted parameter table.
cal$simulator$tmb_model$obj_fn$obj_fn_expr
mp_optimize(cal)
if (suppressPackageStartupMessages(require(broom.mixed))) {
  print(mp_tmb_coef(cal))
}

Description

A ledger is a table with rows that identify specific instances of a functional form used to define a mp_dynamic_model. Ledgers are most commonly created using the mp_join function as in the examples.

Examples

age = mp_index(Age = c("young", "old"))
state = mp_cartesian(
  mp_index(Epi = c("S", "I", "R")),
  age
)
mp_join(
  from = mp_subset(state, Epi = "S"),
  to = mp_subset(state, Epi = "I"),
  by = list(from.to = "Age")
)

Description

An experimental alternative to mp_per_capita_flow that allows users to specify flows using absolute rates instead of per-capita rates.

Usage

mp_absolute_flow(from, to, rate, flow_name = NULL, rate_name = NULL)

Arguments

See Also

mp_per_capita_flow()

Description

Get the adjacency matrix associated with a model specification.

Usage

mp_adjacency(spec, include_inout = FALSE)

Arguments

Value

An adjacency matrix.

Description

Create a one-column ledger (see LedgerDefinition) with rows identifying instances of an aggregation.

Usage

mp_aggregate(index, by = "Group", ledger_column = "group")

Arguments

See Also

Other functions that return ledgers mp_join()

Description

Produce a new index table by taking all possible pairwise combinations of the input tables. This is useful for producing product models that expand model components through stratification.

Usage

mp_cartesian(...)

Arguments

See Also

Other functions that return index tables mp_index(), mp_rename(), mp_subset(), mp_union()

Other functions that take products of index tables and return one index tables mp_linear(), mp_square(), mp_symmetric(), mp_triangle()

Examples

mp_cartesian(
  mp_index(Epi = c("S", "I")),
  mp_index(Age = c("young", "old"))
)

si = mp_index(Epi = c("S", "I"))
age = mp_index(Age = c("young", "old"))
loc = mp_index(City = c("hamilton", "toronto"))
vax = mp_index(Vax = c("unvax", "vax"))
(si
  |> mp_cartesian(age)
  |> mp_cartesian(loc)
  |> mp_cartesian(vax)
)

flow_rates = mp_index(Epi = c("infection", "recovery"))
mp_union(
  mp_cartesian(
    mp_subset(flow_rates, Epi = "infection"),
    age
  ),
  mp_subset(flow_rates, Epi = "recovery")
)

Description

Get a data frame with one row for each change made to each state variable at each time step.

Usage

mp_change_frame(spec)

Arguments

Value

Data frame with two columns: state and change. Each row describes one change.

Examples

("starter_models"
  |> mp_tmb_library("sir", package = "macpan2") 
  |> mp_change_frame()
)

Description

Default Values

Usage

mp_default(model, include_all = FALSE)

mp_default_list(model, include_all = FALSE)

Arguments

Value

A long-format data frame with default values for matrices required as input to model objects. The columns of this output are matrix, row, col, and value. Scalar matrices do not have any entries in the row or col columns.

Functions

• mp_default_list(): List of the default variables as matrices.

Description

This is an 'old' model specification function that was tested out at a workshop. Currently it still drives the engine-agnostic-grammar vignette, but we plan to replace this function with mp_tmb_model_spec and other model specification functions.

Usage

mp_dynamic_model(
  expr_list = ExprList(),
  ledgers = list(),
  init_vecs = list(),
  unstruc_mats = list()
)

Arguments

Description

This is an 'old' function that was tested out at a workshop. Currently it still drives the engine-agnostic-grammar vignette, but we plan to replace this function with mp_simulator.

Usage

mp_dynamic_simulator(
  dynamic_model,
  time_steps = 0L,
  vectors = NULL,
  unstruc_mats = NULL,
  mats_to_save = NULL,
  mats_to_return = NULL,
  params = OptParamsList(0),
  random = OptParamsList(),
  obj_fn = ObjectiveFunction(~0),
  log_file = LogFile(),
  do_pred_sdreport = TRUE,
  tmb_cpp = "macpan2",
  initialize_ad_fun = TRUE,
  ...
)

Arguments

Description

Additional information that can be joined to the output of the tidy.TMB or tidy.stanfit functions in the broom.mixed package.

Usage

mp_effects_descr(model)

mp_add_effects_descr(coef_table, model)

Arguments

Functions

• mp_add_effects_descr(): Convenience function for adding coefficient descriptions from a calibrated model to coef_tables generated by mp_tmb_coef or mp_tmbstan_coef.

Description

Expand a structured model so that it is represented in an unstructured format requiring a more verbose description. Currently, this is only applicable for mp_tmb_model_spec objects that have explicit flows (e.g. mp_per_capita_flow). For such models, mp_expand produces a model with expression lists composed entirely of plain R formulas.

Usage

mp_expand(model)

mp_reduce(model)

Arguments

Functions

• mp_reduce(): Confusingly, mp_reduce and mp_expand are synonyms. Please use mp_expand in new projects, as mp_reduce is available for back-compatibility only.

Examples

sir = mp_tmb_library("starter_models", "sir", package = "macpan2")
print(sir)
print(mp_expand(sir))

Description

Extract the index for a particular dimension in a ledger from a ledger or from an object containing one or more ledgers.

Usage

mp_extract(x, dimension_name)

Arguments

Description

Factor an Index

Usage

mp_factors(index, unpack = c("no", "maybe", "yes"))

Arguments

Description

Return the values of variables after the simulation loop has finished and the final set of expressions have been evaluated.

Usage

mp_final(model)

mp_final_list(model)

Arguments

Functions

• mp_final_list(): Final values formatted as a list of matrices.

Description

Get a data frame where each row represents a flow in a model specification.

Usage

mp_flow_frame(spec, topological_sort = TRUE, loops = "^$")

Arguments

Value

A data frame that gives information provided in calls to mp_per_capita_flow and mp_per_capita_inflow.

Description

A forecaster is an object that can make forecasts. It is constructed from mp_tmb_calibrator, which is an object that can be used to calibrate the parameters of a model – typically by fitting the model to data. The main difference between a calibrator and a forecaster is that the latter runs for more time steps, past the last time step of the calibrator. A forecaster object can be used to generate different types of forecasts using the mp_trajectory family of functions.

Usage

mp_forecaster(
  calibrator,
  forecast_period_time_steps,
  outputs = NULL,
  data = NULL,
  tv = NULL,
  default = list(),
  inits = list()
)

Arguments

Description

Functions Used by an Object for Communicating with a Computational Engine

Usage

mp_functions_used(object)

mp_generates_randomness(object)

Arguments

Value

Character vector of names of functions that are used by object.

Functions

• mp_generates_randomness(): Does the object use functions that generate randomness?

Description

Create a new index with fewer columns to create names for an aggregated vector that is labelled by the input index.

Usage

mp_group(index, by)

Arguments

Description

Make an index table to enumerate model quantity labels by category. These objects generalize and wrap data.frames, where each column is a label category and each row is an index. Indices must contain only letters, numbers, and underscores. Blank empty string entries are allowed, but missing values (NAs) are not.

Usage

mp_index(..., labelling_column_names)

## S3 method for class 'Index'
print(x, ...)

## S3 method for class 'Index'
names(x)

Arguments

Details

Dots are not allowed in indices so that the labels can be inverted to reproduce the original index table (provided that the column names can be retrieved).

It is recommended to use UpperCamelCase for the columns of index tables and single uppercase characters ("S", "I"), all lowercase character strings ("gamma"), and/or snake_case strings ("aging_rate") for indices. This convention helps when reading code that contains references to both column names and indices.

Functions

• print(Index): Print an index.

• names(Index): Get the names of the columns of an index.

See Also

mp_structured_vector()

mp_set_numbers()

Other functions that return index tables mp_cartesian(), mp_rename(), mp_subset(), mp_union()

Examples

state = mp_index(
  Epi = c("S", "I", "S", "I"),
  Age = c("young", "young", "old", "old")
)
print(state)
labels(state)
mp_cartesian(state, mp_index(City = c("hamilton", "toronto")))

# The following index table describes the state variables of the
# model:
sir = mp_index(Epi = c("S", "I", "R"))
print(sir)

# Here, the column `Epi` denotes that the category of these labels is
# epidemiological. There is nothing special about this specific choice of
# category name; we could have also used another name like `Compartment`.

# However, in more complicated models, it is good to think carefully about
# choosing descriptive category names. For example, in an age-structured SIR
# model, we could add an `Age` column to generate an index table as follows:
sir_age = mp_index(
 Epi = rep(c("S", "I", "R"), 2),
 Age = rep(c("young", "old"), each = 3)
)
print(sir_age)

# Here, having the first column in the index table labeled `Compartment` would
# be somewhat misleading, as the compartments aren't actually just "S", "I",
# and "R", they are each of the epidemiological states stratified by the age
# groups "young" and "old".

# This index table could also be generated by first specifying individual index
# tables for the `Epi` and `Age` columns, and then using a `macpan2` product
# function that combines the tables into a single index table:
sir = mp_index(Epi = c("S", "I", "R"))
age = mp_index(Age = c("young", "old"))
prod = mp_cartesian(sir, age)
prod

# The mp_cartesian() function will produce a table with entries that are all
# possible combinations of the individual index tables. The "See Also" section
# of the mp_cartesian() help page catalogues all available product functions.

# We can produce the full labels of model quantities, which are simply
# dot-concatenated indices, one for each entry in the index table, using the
labels(prod)

Description

Return a data frame containing the values of variables at the end of the before phase, right before the simulation loop begins (i.e. right before the during phase).

Usage

mp_initial(model)

mp_initial_list(model)

Arguments

Value

A long-format data frame with initial values for matrices. The columns of this output are matrix, time, row, col, and value. Scalar matrices do not have any entries in the row or col columns. The before phase corresponds to a time value of 0.

Functions

• mp_initial_list(): List of the initial variables as matrices.

Description

Join two or more index tables (see mp_index) to produce a ledger (see LedgerDefinition).

Usage

mp_join(..., by = empty_named_list())

Arguments

Details

When two index tables are passed to ..., mp_join behaves very much like an ordinary inner join. When more than two tables are passed to ..., mp_join iteratively joins pairs of tables to produce a final ledger. For example, if index tables A B, and C are passed to mp_join, an inner join of A and B is performed and the result is joined with C. In each of these successive internal joins. The properties of inner joins ensures that the order of tables does not affect the set of rows in the final table (SW states without proof!).

See Also

Other functions that return ledgers mp_aggregate()

Examples

# When two index tables are passed to `...`, the `by` argument is just a
# character vector of column names on which to join (as in standard R functions
# for joining data frames), or the dot-concatenation of these column names.
# For example,
state = mp_index(
  Epi = c("S", "I", "S", "I"),
  Age = c("young", "young", "old", "old")
)
mp_join(
  from = mp_subset(state, Epi = "S"),
  to = mp_subset(state, Epi = "I"),
  by = "Age"
)
# If there are more than two tables then the `by` argument must be a named
# list of character vectors, each describing how to join the columns of
# a pair of tables in `...`. The names of this list are dot-concatenations
# of the names of pairs of tables in `...`. For example,
rates = mp_index(
  Epi = c("lambda", "lambda"),
  Age = c("young", "old")
)
mp_join(
  from = mp_subset(state, Epi = "S"),
  to = mp_subset(state, Epi = "I"),
  rate = mp_subset(rates, Epi = "lambda"),
  by = list(
    from.to = "Age",
    from.rate = "Age"
  )
)
# If the `by` columns have different names in two tables, then you can
# specify these using formula notation where the left-hand-side
# is a dot-concatenation of columns in the first table and the
# right-hand-side is a dot-concatenation of the columns in the second
# table. For example,
contact = mp_index(
  AgeSusceptible = c("young", "young", "old", "old"),
  AgeInfectious = c("young", "old", "young", "old")
)
mp_join(
  sus = mp_subset(state, Epi = "S"),
  inf = mp_subset(state, Epi = "I"),
  con = contact,
  by = list(
    sus.con = "Age" ~ "AgeSusceptible",
    inf.con = "Age" ~ "AgeInfectious"
  )
)

Description

Convert a scalar binary operator into one that performs a Kronecker product with more convenient row and column names for use with macpan2. The result is numerically identical to base R's %x% and kronecker(), but with cleaner naming of rows and columns.

Usage

mp_kronecker_operator(operator)

Arguments

Value

A Kronecker operator convenient for use with macpan2.

Description

Return a character vector of labels for each row of an index (or a ledger?? FIXME: what does this mean for ledgers??).

Usage

mp_labels(x, labelling_column_names)

Arguments

Description

Create a grid on which to layout the flow diagram of a model specification.

Usage

mp_layout_grid(
  spec,
  east = "",
  south = "^$",
  north = "^$",
  west = "^$",
  loops = north,
  x_gap = 0.3,
  y_gap = 0.3,
  north_south_sep = 0,
  east_west_sep = 0,
  trim_states = character()
)

Arguments

Description

Layout the flow diagram of a model specification so that each row is one of the paths through the model (ignoring loops).

Usage

mp_layout_paths(
  spec,
  sort_paths = TRUE,
  combine_columns = TRUE,
  deduplicate_edges = TRUE,
  loops = "^$",
  ignore = "^$",
  x_gap = 0.3,
  y_gap = 0.3,
  north_south_sep = 0,
  east_west_sep = 0,
  trim_states = character()
)

Arguments

Description

Bundle up several ledgers (see LedgerDefinition) to pass to mp_dynamic_model.

Usage

mp_ledgers(...)

Arguments

Description

TODO: what does this mean?

Usage

mp_linear(x, y_labelling_column_names)

Arguments

See Also

Other functions that take products of index tables and return one index tables mp_cartesian(), mp_square(), mp_symmetric(), mp_triangle()

Description

Lookup a subset or factor index associated with a symbol, and return the index associated with that symbol.

Usage

mp_lookup(index, symbol)

Arguments

Description

Open a browser at the current version of a particular model in an online macpan2 model library.

Usage

mp_model_docs(model_name, macpan_library = "starter_models")

Arguments

Value

This function returns the URL of the library model, but the main purpose is the side-effect of automatically opening a web browser at this URL.

Description

Create a directory with a template model definition.

Usage

mp_model_starter(starter_name, dir)

Arguments

Description

Has an attempt been made to calibrate model parameters through optimization of a likelihood function or posterior density, probably by using mp_optimize?

Usage

mp_opt_attempted(model)

Arguments

Value

Either TRUE or FALSE.

Description

Calibrate a model that has been parameterized, typically by using mp_tmb_calibrator to produce such a model.

Usage

mp_optimize(
  model,
  optimizer = c("nlminb", "optim", "DEoptim", "optimize", "optimise"),
  ...
)

Arguments

Details

The mp_tmb_calibrator models that get passed to mp_optimize will remember both their original default parameters set at the time the model was created, and their currently best parameters that get the updated when mp_optimize is run. The optimization is started from the currently best parameter set. Therefore, if you find yourself in a local minimum you might need to either recreate the model using mp_tmb_calibrator or use optimizer = "DEoptim", which is more robust to objective functions with multiple optima.

Value

The output of the optimizer. The model object is modified and saves the history of optimization outputs. These outputs can be obtained using mp_optimizer_output.

Details on Using Each Type of Optimizer

nlminb

The default optimizer is nlminb. This optimizer uses gradients computed by the Template Model Builder engine of macpan2. This optimizer is efficient at local optimization by exploiting the gradient information computed by automatic differentiation. However, like many nonlinear optimizers, it can struggle if the objective function has multiple optima.

To set control parameters (e.g., maximum number of iterations), you can use the control argument:

mp_optimize(model, "nlminb", control = list(iter.max = 800))

See the nlminb help page for the complete list of control parameters and what the output means.

optim

The optim optimizer lets you choose from a variety of optimization algorithms. The default, method = "Nelder-Mead", does not use second derivatives (compare with the description of nlminb), and so could be less efficient at taking each step. However, we find that it can be somewhat better at getting out of local optima, although the "DEoptim" optimizer is designed for objective functions with multiple optima (see description below).

To set control parameters (e.g., maximum number of iterations), one may use the following.

mp_optimize(model, "nlminb", control = list(maxit = 800))

See the optim help page for the complete list of control parameters and what the output means.

If your model is parameterized by only a single parameter, you'll get a warning asking you to use 'method = "Brent"' or optimizer = 'optimize()'. You can ignore this warning if you are happy with your answer, or can do either of the suggested options as follows.

mp_optimize(model, "optim", method = "Brent", lower = 0, upper = 1.2)
mp_optimize(model, "optimize", interval = c(0, 1.2))

In this case you have to specify lower and upper bounds for the optimization.

DEoptim

The DEoptim optimizer comes from the DEoptim package; you'll need to have that package installed to use this option. It is designed for objective functions with multiple optima. Use this method if you do not believe the fit you get by other methods. The downsides of this method are that it doesn't use gradient information and evaluates the objective function at many different points, so is likely to be much slower than gradient-based optimizers such as the default nlminb optimizer, or optim with method = "BFGS".

Because this optimizer starts from multiple points in the parameter space, you need to specify lower and upper bounds for each parameter in the parameter vector.

mp_optimize(model, "DEoptim", lower = c(0, 0), upper = c(1.2, 1.2))

In this example we have two parameters, and therefore need to specify two values each for lower and upper

optimize and optimise

This optimizer can only be used for models parameterized with a single parameter. You need to specify lower and upper bounds, e.g.

mp_optimize(model, "optimize", c(0, 1.2))

Examples

spec = ("starter_models"
  |> mp_tmb_library("seir", package = "macpan2")
  |> mp_tmb_update(default = list(beta = 0.6))
)
sim = mp_simulator(spec, 50, "infection")

## simulate data to fit to, but remove the start of the
## simulated epidemic in order to make it more difficult
## to fit.
data = mp_trajectory(sim)
data = data[data$time > 24, ]
data$time = data$time - 24

## time scale object that accounts for the true starting time
## of the epidemic (in this example we are not trying to estimate
## the initial number infected, so the starting time strongly
## affects the fitting procedure)
time = mp_sim_offset(24, 0, "steps")

## model calibrator, estimating only the transmission parameter
cal = mp_tmb_calibrator(spec
  , data = data
  , traj = "infection"
  , par = "beta"
  , default = list(beta = 1)
  , time = time
)

## From the starting point at beta = 1 this takes us into a
## local optimum at beta = 0.13, far from the true value of beta = 0.6.
mp_optimize(cal)

## In this case, one-dimensional optimization finds the true value.
mp_optimize(cal, "optimize", c(0, 1.2))

Description

Create a new model specification using parameter values that have been optimized.

Usage

mp_optimized_spec(model, spec_structure = c("original", "modified"))

Arguments

Details

The following is a list of additional model specification structure that was or could have been added to the "modified" specification by mp_tmb_calibrator for the purposes of fitting to data.

• Observed data in a format that can be compared with simulated data.

• Expressions preparing simulated data so that they can be compared with the observed data.

• Distributional parameters for likelihood and prior components of the objective function (e.g., a dispersion parameter for a negatively binomial likelihood component, or a prior standard deviation for a transmission rate parameter).

• Data, parameters, and expressions for modelling time-variation of parameters (e.g., basis functions and weights for a spline determining the time-variation of transmission rate).

Value

A copy of the model specification used to produce the calibrator model, with calibrated default values.

Description

Get the output from an optimizer used in model calibration.

Usage

mp_optimizer_output(model, what = c("latest", "all"))

Arguments

Details

When objects created by mp_tmb_calibrator are successfully passed to mp_optimize, they build up an optimization history. This history is recorded as a list of the output produced by the underlying optimizer (e.g. nlminb). This mp_optimizer_output function returns the latest output by default or the entire history list.

Description

Define the prior distributions for parameters and random effects to be passed to par argument of the mp_tmb_calibrator function.

Usage

mp_par(params, random)

Arguments

Description

Description of Model Parameterization

Usage

mp_parameterization(model, types = c("fixed", "random"))

Arguments

Examples

cal = si_example_object("optimized_calibrator")
mp_parameterization(cal)

Description

Specify different kinds of flows into, out of, and between compartments.

Usage

mp_per_capita_flow(from, to, rate, flow_name = NULL, abs_rate = NULL)

mp_per_capita_inflow(from, to, rate, flow_name = NULL, abs_rate = NULL)

mp_per_capita_outflow(from, rate, flow_name = NULL, abs_rate = NULL)

mp_inflow(to, rate, flow_name = NULL, abs_rate = NULL)

mp_outflow(from, rate, flow_name = NULL, abs_rate = NULL)

Arguments

Details

The examples below can be mixed and matched in mp_tmb_model_spec() to produce compartmental models. The symbols used below must be used in an appropriate context (e.g., if N is used for total population size, then there must be an expression like N ~ S + I + R somewhere in the model or for models with constant population size there must be a default variable, N, with a numerical value).

Functions

• mp_per_capita_inflow(): Only flow into the to compartment, and do not flow out of the from compartment. The from compartment can even be a function of a set of compartments, because it will not be updated. A common construction is mp_per_capita_inflow("N", "S", "birth_rate", "birth") for adding a birth process, which involves the total population size, N, rather than a single compartment.

• mp_per_capita_outflow(): Only flow out of the from compartment, without going anywhere. This is useful for removing individuals from the system (e.g., death). To keep track of the total number of dead individuals one can use mp_per_capita_flow and set to to be a compartment for these individuals (e.g., to = "D").

• mp_inflow(): Only flow into the to compartment. For adding a birth or immigration process.

• mp_outflow(): Only flow out of the from compartment. For adding an absolute removal process that goes to 'nowhere': dangerous! The reason it is dangerous is that this flow can easily lead to negative values of state variables when the rate is high relative to the size of the from compartment. Often mp_per_capita_outflow will be a better choice, given that the size of the outflow will be scaled to the size of the from compartment by measuring rates on a per-capita basis.

See Also

mp_absolute_flow()

Examples

# infection by mass action
# https://github.com/canmod/macpan2/blob/main/inst/starter_models/si
mp_per_capita_flow("S", "I", "beta * I / N", "infection")

# recovery
# https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir
mp_per_capita_flow("I", "R", "gamma", "recovery")

# disease progression with different severity
# https://github.com/canmod/macpan2/blob/main/inst/starter_models/macpan_base
mp_per_capita_flow("E", "I_mild", "alpha * phi"      , "progression_mild")
mp_per_capita_flow("E", "I_sev" , "alpha * (1 - phi)", "progression_sev")

# birth
# https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_demog
mp_per_capita_inflow("N", "S", "nu", "birth")

# death
# https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_demog
mp_per_capita_outflow("S", "mu", "death_S")
mp_per_capita_outflow("I", "mu", "death_I")
mp_per_capita_outflow("R", "mu", "death_R")

# vaccination 
# https://github.com/canmod/macpan2/blob/main/inst/starter_models/shiver
mp_per_capita_flow("S", "V", "(vrate * S/(vrate + S))/S", "vaccination")

# importation
# mp_inflow("I", "delta", "importation")

Description

Return an integer vector of positions of x in table. Currently this is a simple wrapper around match.

Usage

mp_positions(x, table, zero_based = TRUE)

Arguments

Description

Print Objective Function

Usage

mp_print_obj_fn(model)

Arguments

Value

Called to print the objective function for humans to read. Invisibly returns the underlying objective function object that is of limited utility.

Description

Print Model Specification

Usage

mp_print_spec(model)

mp_print_before(model)

mp_print_during(model)

mp_print_after(model)

Arguments

Functions

• mp_print_before(): Print just the expressions executed before the simulation loop.

• mp_print_during(): Print just the expressions executed during each iteration of the simulation loop.

• mp_print_after(): Print just the expressions executed after the simulation loop.

Description

Pass the output of this function to the tv argument of mp_tmb_calibrator to model time variation of a parameter with flexible radial basis functions.

Usage

mp_rbf(
  tv,
  dimension,
  initial_weights,
  seed,
  prior_sd = 1,
  fit_prior_sd = TRUE,
  sparse_tol = 0.01
)

Arguments

Description

Uses readRDS to read in a saved model specification created using a function like mp_tmb_model_spec, and updates this specification using mp_version_update so that it is compatible with the installed version of macpan2. To save a model specification, just use the base R function saveRDS.

Usage

mp_read_rds(filename)

Arguments

Description

Extract the index used as a reference for generating position vectors.

Usage

mp_reference(x, dimension_name)

Arguments

Description

Rename Index Columns

Usage

mp_rename(x, ...)

Arguments

See Also

Other functions that return index tables mp_cartesian(), mp_index(), mp_subset(), mp_union()

Description

Collects information from the headers of the README files in model directories and returns the results as a data frame

Usage

mp_show_models(
  dir = system.file("starter_models", package = "macpan2"),
  show_missing = FALSE,
  for_markdown = FALSE
)

show_models(
  dir = system.file("starter_models", package = "macpan2"),
  show_missing = FALSE,
  for_markdown = FALSE
)

mp_list_models(dir = system.file("starter_models", package = "macpan2"))

Arguments

Value

a data frame containing entries Directory (model directory), Title (model title), Description (short description)

Functions

• show_models(): Synonym for mp_show_models, which is preferred. Present for back-compatibility.

• mp_list_models(): Return a character vector containing model names, instead of a data frame with more information about each model.

Examples

mp_show_models(show_missing = TRUE)

Description

Set the simulation bounds (start time and end time) for a calibration. This is used to override the default simulation bounds taken from the observed data passed to mp_tmb_calibrator.

Usage

mp_sim_bounds(sim_start, sim_end, time_scale = "steps", time_column = "time")

Arguments

Value

An object to be passed to the time argument of mp_tmb_calibrator.

See Also

mp_sim_offset()

Description

Offset the starting and ending times of the simulation, from the start and end time of the data used in calibration. This is used to override the default offsets of zero taken from the observed data passed to mp_tmb_calibrator.

Usage

mp_sim_offset(
  sim_start_offset,
  sim_end_offset,
  time_scale = "steps",
  time_column = "time"
)

Arguments

Value

An object to be passed to the time argument of mp_tmb_calibrator.

See Also

mp_sim_bounds()

Description

Construct a simulator from a model specification object.

Usage

mp_simulator(model, time_steps, outputs, default = list(), inits = list())

Arguments

Description

Slice an index

Usage

mp_slices(index, unpack = c("no", "maybe", "yes"))

Arguments

Description

Self Cartesian Product

Usage

mp_square(x, suffixes = c("A", "B"))

Arguments

See Also

Other functions that take products of index tables and return one index tables mp_cartesian(), mp_linear(), mp_symmetric(), mp_triangle()

Description

Data frame giving states that per-capita flow rates directly depend on. This is intended for plotting diagrams and not for mathematical analysis, in that it does not describe indirect dependence for flow rates on state variables.

Usage

mp_state_dependence_frame(spec)

Arguments

Description

This documentation was originally in mp_index() and should be cleaned up See issue #131. Also this is an experimental feature.

Usage

mp_structured_vector(x, ...)

mp_set_numbers(vector, ...)

Arguments

Functions

• mp_set_numbers(): Update numerical values of a structured vector. TODO: details on syntax.

Examples

state = mp_index(
  Epi = c("S", "I", "S", "I"),
  Age = c("young", "young", "old", "old")
)
state_vector = (state
  |> mp_structured_vector()
  |> mp_set_numbers(Epi = c(S = 1000))
  |> mp_set_numbers(Epi = c(I = 1), Age = "old")
)
print(state_vector)

Description

Take a subset of the rows of an index table (see mp_index) to produce another index table. The mp_subset function gives rows that match a certain criterion and mp_setdiff gives rows that do not match.

Usage

mp_subset(x, ...)

mp_setdiff(x, ...)

Arguments

See Also

Other functions that return index tables mp_cartesian(), mp_index(), mp_rename(), mp_union()

Description

Symmetric Self Cartesian Product

Usage

mp_symmetric(x, y_labelling_column_names, exclude_diag = TRUE)

Arguments

See Also

Other functions that take products of index tables and return one index tables mp_cartesian(), mp_linear(), mp_square(), mp_triangle()

Description

Time Scale

Usage

mp_time_scale(start, end, time_step_scale = c("steps", "daily", "weekly"), ...)

Arguments

Description

Get the result of TMB::MakeADFun underlying a TMB-based model in macpan2.

Usage

mp_tmb(model)

Arguments

Description

Construct an object that can get used to calibrate an object produced by mp_tmb_model_spec or mp_tmb_library, and possibly modified by mp_tmb_insert or mp_tmb_update. Note that the output of this function is not a model that has been calibrated to data, but rather an object that contains a model that can be calibrated. To actually attempt a calibration one needs the mp_optimize function (see examples below).

Usage

mp_tmb_calibrator(
  spec,
  data = empty_trajectory,
  traj = character(),
  par = character(),
  tv = character(),
  outputs = traj,
  default = list(),
  inits = list(),
  time = NULL,
  save_data = TRUE,
  optimize = FALSE
)

Arguments

Examples

spec = mp_tmb_library("starter_models", "sir", package = "macpan2")
sim = mp_simulator(spec, 50, "infection")
data = mp_trajectory(sim)
cal = mp_tmb_calibrator(
    spec
  , data
  , traj = "infection"
  , par = "beta"
  , default = list(beta = 0.25)
)
mp_optimize(cal)
if (suppressPackageStartupMessages(require(broom.mixed))) {
  print(mp_tmb_coef(cal))
}

Description

TMB Model Coefficient Table

Usage

mp_tmb_coef(model, back_transform = TRUE, ...)

Arguments

Value

A data frame that describes the fitted coefficients.

Description

Create a list of expressions for defining a compartmental model in TMB.

Usage

mp_tmb_expr_list(
  before = list(),
  during = list(),
  after = list(),
  .simulate_exprs = character(0L)
)

Arguments

Value

Object of class ExprList with the following methods.

Methods

• ⁠$data_arg(...)⁠: Return the following components of the data structure to pass to C++.

  • expr_output_id – Indices into the list of matrices identifying the matrix being produced.

  • expr_sim_block – Identified whether or not the expression should be evaluated inside a simulate macro within TMB.

  • expr_num_p_table_rows – Number of rows associated with each expression in the parse table (⁠p_table_*⁠)

  • eval_schedule – Vector giving the number of expressions to evaluate in each phase (before, during, or after) of the simulation.

  • p_table_x – Parse table column giving an index for looking up either function, matrix, or literal.

  • p_table_n – Parse table column giving the number of arguments in functions.

  • p_table_i – Parse table column giving indices for looking up the rows in the parse table corresponding with the first argument of the function.

Method Arguments

• ...: Character vector containing the names of the matrices in the model.

Description

Covariance of Fixed Effect Estimates

Usage

mp_tmb_fixef_cov(model)

Arguments

Value

A covariance matrix.

Description

Insert, update, or delete elements of a TMB model spec, produced using mp_tmb_library or mp_tmb_model_spec. The only difference between mp_tmb_insert and mp_tmb_update is that the former shifts the positions of existing expressions to make room for the new expressions, whereas the latter overwrites existing expressions using the new expressions. The treatment of new default values and integers is the same. The examples below clarify this difference. Note that mp_tmb_delete does not contain an expressions argument, because it is not necessary to specify new expressions in the case of deletion.

Usage

mp_tmb_insert(
  model,
  phase = "during",
  at = 1L,
  expressions = list(),
  default = list(),
  inits = list(),
  integers = list(),
  must_save = character(),
  must_not_save = character(),
  sim_exprs = character()
)

mp_tmb_update(
  model,
  phase = "during",
  at = 1L,
  expressions = list(),
  default = list(),
  inits = list(),
  integers = list(),
  must_save = character(),
  must_not_save = character(),
  sim_exprs = character()
)

mp_tmb_delete(
  model,
  phase,
  at,
  default = character(),
  integers = character(),
  must_save = character(),
  must_not_save = character(),
  sim_exprs = character()
)

Arguments

Details

These modifications do not update the model specification in-place. Instead the output of mp_tmb_insert, mp_tmb_update, and mp_tmb_delete define a new model specification and should be saved if you want to use the new model (e.g., new_model = mp_tmb_insert(model, ...)).

Value

A new model spec object with updated and/or inserted information.

Examples

si = mp_tmb_library("starter_models", "si", package = "macpan2")
print(si)

## Update the mixing process to include 
## optional phenomenological heterogeneity.
## We need mp_tmb_update here so that 
## the previous infection expression is
## overwritten.
mp_tmb_update(si, phase = "during"
  , at = 1
  , expressions = list(infection ~ beta * I * (S/N)^zeta)
  , default = list(zeta = 1)
)

## Parameterize with log_beta in place of beta.
## We need mp_tmb_insert here so that the
## existing expression for computing the initial
## number of susceptible indiviudals is not
## overwritten.
mp_tmb_insert(si, phase = "before"
  , at = 1
  , expressions = list(beta ~ exp(log_beta))
  , default = list(log_beta = log(0.5))
)

Description

Insert Back Transformations of Model Parameters

Usage

mp_tmb_insert_backtrans(
  model,
  variables = character(),
  transformation = mp_log
)

mp_tmb_implicit_backtrans(model, variables = character())

Arguments

Value

A new model specification object with expressions for the untransformed (or back transformed) parameters at the beginning of the before phase. The transformed version of the parameter is also added to the defaults and are identified with a prefixed name (e.g., log_beta if beta is log transformed).

Functions

• mp_tmb_implicit_backtrans(): Insert parameter transformations implicitly by pre-pending the name of the transformations in front of the names of the variables that you ask for (e.g., "log_case_reports") even if this variable is not in the model as long as the base name (e.g., "case_reports") is.

See Also

mp_tmb_insert_trans()

Examples

init_si = ("starter_models"
  |> mp_tmb_library("si", package = "macpan2")
  |> mp_tmb_insert_backtrans("beta", mp_log)
  |> mp_initial_list()
)
print(init_si$log_beta)
print(log(init_si$beta))

Description

Insert GLM Time Variation

Usage

mp_tmb_insert_glm_timevar(
  model,
  parameter_name,
  design_matrix,
  timevar_coef,
  link_function = mp_log,
  matrix_coef_name = sprintf("matrix_coef_%s", parameter_name),
  matrix_row_name = sprintf("matrix_row_%s", parameter_name),
  matrix_col_name = sprintf("matrix_col_%s", parameter_name),
  linear_pred_name = sprintf("linear_pred_%s", parameter_name),
  timeseries_name = sprintf("timeseries_%s", parameter_name),
  timevar_coef_name = sprintf("time_var_%s", parameter_name),
  time_index_name = sprintf("time_index_%s", parameter_name),
  sparsity_tolerance = 0,
  engine_function = c("sparse_mat_mult", "group_sums")
)

Arguments

Description

Insert Log Linear Model of Time Variation (Experimental)

Usage

mp_tmb_insert_log_linear(
  model,
  parameter_name,
  design_matrices,
  time_var_parameters,
  window_names = names(time_var_parameters),
  baseline_functions = c(list(TimeVarBaselineParameter()),
    rep(list(TimeVarBaselineNumeric(0)), length(design_matrices) - 1)),
  link_functions = rep(list(mp_identity), length(design_matrices)),
  full_series_name = sprintf("time_var_output_%s", parameter_name),
  baseline_names = sprintf("baseline_%s", window_names),
  matrix_coef_names = sprintf("matrix_coef_%s", window_names),
  matrix_row_names = sprintf("matrix_row_%s", window_names),
  matrix_col_names = sprintf("matrix_col_%s", window_names),
  linear_pred_names = sprintf("linear_pred_%s", window_names),
  time_var_names = sprintf("time_var_%s", window_names),
  time_index_name = sprintf("time_index_%s", parameter_name),
  sparsity_tolerance = 0
)

Arguments

Description

A version of mp_tmb_insert making it more convenient to transform an incidence variable into a reports variable, which accounts for reporting delays and under-reporting. This new reports variable is a convolution of the simulation history of an incidence variable with a kernel that is proportional to a Gamma distribution of reporting delay times.

Usage

mp_tmb_insert_reports(
  model,
  incidence_name,
  report_prob,
  mean_delay,
  cv_delay,
  reports_name = sprintf("reported_%s", incidence_name),
  report_prob_name = sprintf("%s_report_prob", incidence_name),
  mean_delay_name = sprintf("%s_mean_delay", incidence_name),
  cv_delay_name = sprintf("%s_cv_delay", incidence_name)
)

Arguments

Description

Insert Basic Transformations of Model Variables

Usage

mp_tmb_insert_trans(model, variables = character(), transformation = mp_log)

mp_tmb_implicit_trans(model, variables = character())

Arguments

Value

A new model specification object with expressions for the transformed variables at the end of the simulation loop. The transformed variables are identified with a prefixed name (e.g., log_incidence if incidence is log transformed).

Functions

• mp_tmb_implicit_trans(): Insert variable transformations implicitly by pre-pending the name of the transformations in front of the names of the variables that you ask for (e.g., "log_case_reports") even if this variable is not in the model as long as the base name (e.g., "case_reports") is.

See Also

mp_tmb_insert_backtrans()

Examples

("starter_models"
  |> mp_tmb_library("si", package = "macpan2")
  |> mp_tmb_insert_trans("infection", mp_log)
  |> mp_simulator(time_steps = 5L, outputs = "log_infection")
  |> mp_trajectory()
)

Description

Get a TMB model specification from a model library.

Usage

mp_tmb_library(..., package = NULL, alternative_specs = FALSE)

mp_tmb_entire_library()

mp_official_library(model_name)

Arguments

Details

This function executes the R code in the tmb.R file of a directory in a model library. To be a valid model, this file should produce an object called spec (containing a model specification produced by the mp_tmb_model_spec() function) or specs (containing a list of such specifications). This mp_tmb_library() function returns this spec or specs object (see the alternative_specs argument), but does not expose any other object produced by tmb.R.

Functions

• mp_tmb_entire_library(): List of model specifications; one for each model in the library.

• mp_official_library(): Get a model specification from the official macpan2 starter_models library.

See Also

mp_show_models()

Examples

mp_tmb_library(
    "starter_models"
  , "si"
  , package = "macpan2"
)

Description

Specify a simulation model in the TMB engine. A detailed explanation of this function is covered in vignette("quickstart").

Usage

mp_tmb_model_spec(
  before = list(),
  during = list(),
  after = list(),
  default = list(),
  inits = list(),
  integers = list(),
  must_save = character(),
  must_not_save = character(),
  sim_exprs = character(),
  state_update = c("euler", "rk4", "discrete_stoch", "hazard", "rk4_old",
    "euler_multinomial")
)

Arguments

Details

Expressions in the before, during, and after lists can be standard R formula objects for defining variables in the model. These formulas must have a left hand side that gives the name of the (possibly matrix-valued) variable being updated, and a right hand side giving an expression containing only (1) the names of quantities in the model, (2) numerical literals (e.g., 3.14), or (3) functions defined in the TMB engine (described in engine_functions). For example, the expression N ~ S + I + R updates the value of N to be the sum of the variables S, I, and R.

Names can be provided for the components of the before, during, and after lists, and these names do not have to be unique. These names are used by the sim_exprs argument.

Examples

## A simple SI model.
spec = mp_tmb_model_spec(
    during = mp_per_capita_flow("S", "I", "beta * I / N", "infection")
  , default = list(N = 100, S = 99, I = 1, beta = 0.2)
)
(spec 
  |> mp_simulator(time_steps = 5L, output = "infection") 
  |> mp_trajectory()
)

## The `~` can be used for flexibly defining dynamical variables.
spec2 = mp_tmb_model_spec(
    during = list(
          force_of_infection ~ beta * I / N
        , mp_per_capita_flow("S", "I", "force_of_infection", "infection")
    )
  , default = list(N = 100, S = 99, I = 1, beta = 0.2)
)
(spec2
  |> mp_simulator(time_steps = 5L, output = "force_of_infection") 
  |> mp_trajectory()
)

## The `before` argument can be used to pre-compute quantities before
## the simulation loop begins. Here we compute `S` from `N` and `I`,
## instead of specifying `S` in the `default` list. The potential
## benefit here is that one could make `I` a parameter to be fitted,
## while ensuring consistent values for `S`.
spec3 = mp_tmb_model_spec(
    before = S ~ N - I
  , during = mp_per_capita_flow("S", "I", "beta * I / N", "infection")
  , default = list(N = 100, I = 1, beta = 0.2)
)
(spec3 
  |> mp_simulator(time_steps = 5L, output = "infection") 
  |> mp_trajectory()
)

Description

Value of the Objective Function of a Model

Usage

mp_tmb_objective(
  model,
  parameter_updates = list(),
  baseline = c("recommended", "default", "optimized")
)

Arguments

Description

Use TMB::tmbprofile to compute the profile likelihood of a calibrator produced using mp_tmb_calibrator.

Usage

mp_tmb_profile(model, param, ...)

Arguments

Value

The output of TMB::tmbprofile.

Description

Leverages the tmbstan and broom.mixed packages to generate MCMC-based coefficient tables.

Usage

mp_tmbstan_coef(model, tmbstan_args = list(), ...)

Arguments

Description

Specify a set of trajectories to fit. The output of this function is intended to be passed to the traj argument of mp_tmb_calibrator.

Usage

mp_traj(likelihood = empty_named_list(), condensation = empty_named_list())

Arguments

Description

Return simulations of the trajectory of the output variables of a dynamical model simulator. To see this functionality in context, please see vignette("quickstart").

Usage

mp_trajectory(model, include_initial = FALSE)

mp_trajectory_par(
  model,
  parameter_updates = list(),
  include_initial = FALSE,
  include_final = FALSE,
  baseline = c("recommended", "default", "optimized")
)

mp_trajectory_sd(
  model,
  conf.int = FALSE,
  conf.level = 0.95,
  include_initial = FALSE,
  back_transform = TRUE
)

mp_trajectory_ensemble(model, n, probs = c(0.025, 0.975))

mp_trajectory_sim(model, n, probs = c(0.025, 0.25, 0.5, 0.75, 0.975))

mp_trajectory_replicate(
  model,
  n,
  parameter_updates = list(),
  include_initial = FALSE,
  include_final = FALSE,
  baseline = c("recommended", "default", "optimized")
)