Defining NLME models in Pumas

We provide two model types, PumasModel and PumasEMModel. The PumasModel type provides great flexibility, allowing a wide variety of random effects distributions and error models and the best support for different approximation methods, excluding only SAEM. The PumasEMModel type is more structured, supporting just a subset of the models expressible by a PumasModel, but uses this structure to enable the optimizations required to make SAEM efficient.

Instances of PumasModel and PumasEMModel types are easily constructed using a macro based domain-specific language (DSL) using the macros @model and @emmodel respectively.

The @model macro interface

To define an NLME model in Pumas, a PumasModel, we use the @model macro. We can define the simplest model of them all, the empty model, as follows:

@model begin

end

This creates a model with no parameters, no covariates, no dynamics, ..., nothing! To populate the model, we need to include some other possible blocks inside the @model begin ... end block. The available blocks are:

• @param: fixed effects and distributional parameters specifications
• @random: random effects specifications
• @covariates: covariate names
• @pre: pre-processing variables for the dynamic system and statistical specification
• @dosecontrol: specification of any dose control parameters present in the model
• @vars: shorthand notation
• @init: initial conditions for the dynamic system
• @dynamics: dynamics of the model (or a @reactions block can be used instead)
• @derived: statistical modeling of dependent variables
• @observed: model information to be stored in the model solution
• @options: specifies model options, including checklinear, inplace, and subject_t0.

In general, it is common to specify blocks in this order (first @param, then @random etc...), where often definitions in some blocks are used by other blocks further down the list.

@param: Population parameters

The population parameters are specified in the @param block. In a PumasModel, variables that enter the model can either be defined in terms of the domain they come from or their distribution if they're random variables. Variables defined by their domain are specified by an in (or ∈, written via \in<TAB>) statement that connects a parameter name and a domain. Random variables are specified by a ~ statement that connects a name with a distribution. For example, to specify θ as a real scalar in a model, one would write:

@model begin
  @param begin
    θ ∈ RealDomain(; lower = 0.0, upper = 17.0)
  end
end

# output

PumasModel
  Parameters: θ
  Random effects:
  Covariates:
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived:
  Observed:

which creates a model with a parameter that has a lower and upper bound on the allowed values.

Different domains are available for different purposes. Their names and purposes are:

• RealDomain for scalar parameters
• VectorDomain for vectors
• PDiagDomain for positive definite matrices with diagonal structure
• PSDDomain for general positive semi-definite matrices

Different domains can be used when we want to have our parameters be scalars or vectors (RealDomain vs VectorDomain) or have certain properties (PDiagDomain and PSDDomain). The simplest way of specifying a model is in terms of all scalar parameters:

@model begin
  @param begin
    θCL ∈ RealDomain(; lower = 0, upper = 50)
    θV ∈ RealDomain(; lower = 0, upper = 500)
    ω²η ∈ RealDomain(; lower = 0, upper = 20)
  end
end

where we defined three separate scalar variables. The parameter θCL is the population clearance, the parameter θV is the volume of distribution, and the parameter ω²η is the variance of a scalar (univariate) random effect. Note that to input the last variable use\omega<TAB>\^2<TAB>\eta<TAB>.

An equivalent model could be instead specified using vector and matrix type domains using something like the following:

@model begin
  @param begin
    θ ∈ VectorDomain(2; lower = [0, 0], upper = [50, 500])
    Ωη ∈ PDiagDomain(1) # no lower or upper keywords!
  end
end

Notice, that we collapsed the two parameters θCL and θV into a single vector θ, and if we want to use the elements in the model you will have to use indexing θ[1] for θCL and θ[2] for θV. It is also necessary to specify the dimension of the vector which is two in this case. The PDiagDomain domain type is special. It makes Ωη have the interpretation of a matrix type, specifically a diagonal matrix. Additionally, it tells Pumas that when fiting, the multivariate parameter (a single scalar in this case) should be kept positive definite. The obvious use case here is variance-covariance matrices, and specifically it is useful for random effect vectors where each random effect is independent of the other. We will get back to this below.

Finally, we have the PSDDomain. This is different from PDiagDomain mainly by representing a "full" variance-covariance matrix. This domain allows the model to have correlation between random effects.

@model begin
  @param begin
    θ ∈ VectorDomain(2; lower = [0, 0], upper = [50, 500])
    Ωη ∈ PSDDomain(1) # no lower upper!
  end
end

As stated above, in addition to the non-random domains, we may specify parameters via a distribution. This allows us to specify parameter in terms of their priors. For example, a model with a parameter that has a multivariate normal (MvNormal) prior can be defined as:

μ_prior = [0.1, 0.3]
Σ_prior = [
  1.0 0.1
  0.1 3.0
] #Observe that there are no commas (,) when constructing a matrix in Julia

@model begin
  @param begin
    θ ~ MvNormal(μ_prior, Σ_prior)
  end
end

A prior can be wrapped in a Constrained(prior; lower=lv, upper=uv) to constrain values to be between lv and uv. See ?Constrained for more details.

@random: Random effects

A key strength of the NLME approach for pharmacometrics is its ability to handle random effects for subject variability and other factors. We just saw that @param was used to specify fixed effects. The appropriate block for specifying random effects is simply called @random. In a PumasModel, the parameters specified in this block can be scalars or vectors just as fixed parameters, but they are always defined by the distribution they are assumed to follow.

Similarly to parameters with priors illustrated above, random effect parameters of a PumasModel are defined by a ~ (read: distributed as) expression. These are used inside a @random block as follows:

@model begin
  @param begin
    ωη ∈ RealDomain(; lower = 0)
  end
  @random begin
    η ~ Normal(0, ωη)
  end
end

# output

PumasModel
  Parameters: ωη
  Random effects: η
  Covariates:
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived:
  Observed:

Here, we defined a variability parameter ωη to parameterize the standard deviation of the univariate Normal distribution of the η. We put a lower bound of 0.0001 because a standard deviation cannot be negative, and a standard deviation of exactly zero would lead to a degenerate distribution. We always advise putting bounds on variables whenever possible. If the parameter was a variance as in a previous example, we advise using the Unicode ² (\^2 + Tab) as part of the name, though this is optional. The Normal distribution Distributions requires two positional arguments: the mean (here: 0) and the standard deviation, here ωη. For more details type ?Normal in the REPL.

Non-Gaussian random effects can also be defined as such:

@model begin
  @param begin
    α ∈ RealDomain(; lower = 0)
    β ∈ RealDomain(; lower = 0)
  end
  @random begin
    η ~ Beta(α, β)
  end
end

# output

PumasModel
  Parameters: α, β
  Random effects: η
  Covariates:
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived:
  Observed:

The following are all the continuous univariate distributions available for use in a PumasModel's @random block. For more details on any distribution, type ? in the REPL followed by the distribution name, e.g. ?Beta.

• Arcsine
• Beta
• BetaPrime
• Biweight
• Cauchy
• Chi
• Chisq
• Cosine
• Epanechnikov
• Erlang
• Exponential
• FDist
• Frechet
• Gamma
• GeneralizedExtremeValue
• GeneralizedPareto
• Gumbel
• InverseGamma
• InverseGaussian
• Kolmogorov
• LocationScale
• Logistic
• LogitNormal
• LogNormal
• NormalCanon
• PGeneralizedGaussian
• Rayleigh
• Semicircle
• SymTriangularDist
• TDist
• Triweight
• Uniform
• Weibull

It is, of course, possible to have as many univariate (or multivariate) random effects as you want. For example:

@model begin
  @param begin
    ω²η ∈ VectorDomain(3; lower = 0)
  end
  @random begin
    η1 ~ Normal(0, sqrt(ω²η[1]))
    η2 ~ Normal(0, sqrt(ω²η[2]))
    η3 ~ Normal(0, sqrt(ω²η[3]))
  end
end

Notice the use of indexing into the ω²η parameter that is now a vector of variances. Other ways of parameterizing random effects include vector (multivariate) distributions:

@model begin
  @param begin
    ω²η ∈ VectorDomain(3, lower = 0)
  end
  @random begin
    η ~ MvNormal(sqrt.(ω²η))
  end
end

where η will have a diagonal variance-covariance structure because we input a vector of standard deviations to MvNormal. This can also be achieved using the PDiagDomain as we saw earlier, and then we don't have to worry about the lower keyword:

@model begin
  @param begin
    Ωη ∈ PDiagDomain(3)
  end
  @random begin
    η ~ MvNormal(Ωη)
  end
end

Do notice that Ωη is not to be considered a vector here, but an actual diagonal matrix. Hence, Ωη is now the (diagonal) variance-covariance matrix, not a vector of standard deviations. The @random block is the same if we allow full covariance structure:

@model begin
  @param begin
    Ωη ∈ PSDDomain(3)
  end
  @random begin
    η ~ MvNormal(Ωη)
  end
end

For cases where you have several random effects with the exact same distribution, such as between-occasion-variability (BOV), it is convenient to construct a single vector η that has diagonal variance-covariance structure with identical variances down the diagonal. This can be achieved using a special MvNormal constructor that takes in the dimension and the standard deviation:

@model begin
  @param begin
    ω²η ∈ RealDomain(; lower = 0)
  end
  @random begin
    η ~ MvNormal(4, sqrt(ω²η))
  end
end

You could use four scalar η's as shown above, but for BOV it is useful to encode the occasions using integers 1, 2, 3, ..., N and simply index into η using η[OCC] where OCC is the occasion covariate.

@covariates: covariate names

The covariates in the remaining model blocks have to be specified in the @covariates block. This information is used to generate efficient code for expanding covariate information from each subject when fitting the model or evaluating likelihood contributions from observations. The format is simply to either use a one-liner:

@model begin
  @covariates weight age OCC
end

# output

PumasModel
  Parameters:
  Random effects:
  Covariates: weight, age, OCC
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived:
  Observed:

or as a block:

@model begin
  @covariates begin
    weight
    age
    OCC
  end
end

# output

PumasModel
  Parameters:
  Random effects:
  Covariates: weight, age, OCC
  Dynamical system variables:
  Dynamical system type: No dynamical model
  Derived:
  Observed:

The block form is mostly useful if there are a lot of covariates. Otherwise, the one-liner is preferred.

@pre: Pre-processing of input to dynamics and derived

Before we move to the actual dynamics of the model (if there are any) and the statistical model of the observed variables we need to do some preprocessing of parameters and covariates to get our rates and variables ready for our ODEs or distributions. This is done in the @pre block.

In the @pre block all calculations are written as if they happen at some arbitrary point in time t. Let us see an example:

@model begin
  # Fixed parameters
  @param begin
    θCL ∈ RealDomain(; lower = 0, upper = 20)
    θV ∈ RealDomain(; lower = 0, upper = 91)
    θbioav ∈ RealDomain(; lower = 0, upper = 1)
    ω²η ∈ RealDomain(; lower = 0)
  end

  # Random parameters
  @random begin
    η ~ MvNormal(4, sqrt(ω²η))
  end

  # Covariate enumeration
  @covariates weight age OCC

  # Preprocessing of input to dynamics and derived
  @pre begin
    CL = θCL * sqrt(weight) / age + η[OCC]
    V = θV
  end
end

We see that when we assign the right-hand side to CL, it involves θCL, weight, age, and the occasion counter, OCC. These might all be recorded as time-varying, especially the last one. The first line of @pre then means that whenever CL is referenced in the dynamic model or in the statistical model it will have been calculated with the covariates evaluated at the appropriate time.

If the dynamic parameters such as clearance or volume are modeled as continuous functions of time (as opposed to piece-wise constant in the covariate case), it is possible to use the reserved named t inside @pre. An example of this kind of model would be a model with enterohepatic recirculation

pk = @model  begin
   @param begin
     tvcl        ∈ RealDomain(lower=0)
     tvvc        ∈ RealDomain(lower=0)
     tvvp        ∈ RealDomain(lower=0)
     tvQ         ∈ RealDomain(lower=0)
     tvka        ∈ RealDomain(lower=0)
     tvklg       ∈ RealDomain(lower=0)
     tvτ         ∈ RealDomain(lower=0)
     Ω           ∈ PDiagDomain(7)
     σ²_prop     ∈ RealDomain(lower=0)
   end

   @random begin
     η           ~ MvNormal(Ω)
   end

   @pre begin
     Cl          = tvcl * exp(η[1])
     Vc          = tvvc * exp(η[2])
     Vp          = tvvp * exp(η[3])
     Q           = tvQ * exp(η[4])
     Ka          = tvka * exp(η[5])
     Klg         = tvklg * exp(η[6])
     τ           = tvτ * exp(η[7])
     Kempt       = (t>10 && t<(10+τ))*(1/τ)
   end
end

where Kempt is the so-called bile emptying rate constant that directly depends on time.

Only variables explicitly defined in @pre can be used in the @dynamics block below. This means that even parameters that require no further pre-processing before they're using in the model will have to be assigned a name in @pre. For example, the following example is the appropriate way to specify a parameter V that is going to be used in the model:

@model begin
  @param begin
    θV ∈ RealDomain(; lower = 0, upper = 91)
  end

  @pre begin
    V = θV
  end
end

If you omit the definition of V and try to use θV directly in the model, you will get an error!

@dosecontrol: Dose Control Parameters

@dosecontrol begin
  bioav = (; Depot = θbioav, Central = 0.4)
end

The @dosecontrol block controls the dose control parameters (DCPs) of the model. The line sets bioavailability for a Depot compartment to a parameter in our model θbioav, and bioavailability of another compartment Central to a fixed value of 0.4. If a compartment is not mentioned in the NamedTuple it will be set to 1.0. Other DCPs are: lags, rate, and duration. For more information on these parameters, see the Dose Control Parameters (DCP) page.

The @dosecontrol block does not share computations with @pre. The available variables are the same as in @pre. These are those defined in @param, @random, @covariates, the current time t, and events of the subject.

@vars: Short-hand notation

Suppose we have a model with a dynamic variable Central and a volume of distribution V. You can define shorthand notation for the implied plasma concentration to be used elsewhere in the model in @vars:

@model begin
  ...
  @vars begin
    conc = Central / V
  end
end

While some users find @vars useful we advise users to use it with caution. Shorthand notation involving dynamic variables might make the @dynamics block harder to read. Shorthand notation that doesn't involve dynamic variables should rather just be specified in @pre. @vars is specific to PumasModels.

@init: Initializing the dynamic system

This block defines the initial conditions of the dynamical model in terms of the parameters, random effects, and pre-processed variables. It is defined by a series of equality (=) statements. For example, to set the initial condition of the response dynamical variable to be the value of the 2nd term of the parameter θ, we would use the syntax:

@model begin
  @param begin
    θ ∈ VectorDomain(3; lower = [0, 0, 1], upper = [3, 1, 4])
  end
  @init begin
    Depot = θ[2]
  end
end

Any variable omitted from this block is given the default initial condition of 0. If the block is omitted, then all dynamical variables are initialized at 0.

Note that the special character combination := can be used to define intermediate statements that will not be carried outside the block. This means that all the resulting data workflows from this model will not contain the intermediate variables defined with :=. For example the result from an inspect call would not contain the intermediate variables that were defined with :=.

@dynamics: The dynamic model

The @dynamics block defines the nonlinear function from the parameters to the derived variables via a dynamical (differential equation) model. It can be specified either by an analytical solution type, an ordinary differential equation (ODE), or a combination of the two.

The analytical solutions are defined in the Analytical Solutions and Differential Equations page and can be invoked via the name. For example:

@model begin
  @param begin
    θCL ∈ RealDomain(; lower = 0, upper = 10)
    θVc ∈ RealDomain(; lower = 0, upper = 10)
  end

  @pre begin
    CL = θCL
    Vc = θVc
  end

  @dynamics Central1
end

defines the dynamical model as the one compartment model represented by Central1. The model has two required parameters: CL and Vc. These have to be defined in @pre when this model is used. All models with analytical solutions have the required parameters listed in their docstring which can be seen by typing ?Central1 in the REPL. Alternatively, it is listed in the documentation on the Analytical Solutions page.

For a system of ODEs that has to be numerically solved, the dynamical variables are defined by their derivative expression. A derivative expression is given by a variable's derivative (specified by ') and an equality (=). For example, the following defines a model equivalent to the model above:

@model begin
  @param begin
    θCL ∈ RealDomain(; lower = 0, upper = 10)
    θVc ∈ RealDomain(; lower = 0, upper = 10)
  end

  @pre begin
    CL = θCL
    Vc = θVc
  end

  @dynamics begin
    Central' = -CL / Vc * Central
  end
end

Variable aliases defined in the @vars are accessible in this block. Additionally, the variable t is reserved for the solver time if you want to use something like sin(t) in your model formulation.

Note that any Julia function defined outside the @model block can be invoked in the @dynamics block.

@reactions: The dynamics model in chemical arrow notation

The @reactions block is a stand-in replacement for the @dynamics block, and the two serve the same purpose. However, in the @reactions block, the modeled interactions are represented through a chemical arrow notation. For example, (k_f, k_b), A + B <--> C represents the association and dissociation between A and B to form the complex C at forward rate k_f and backward rate k_b. The equivalent ODE representation would be:

\[\begin{align*} \frac{dA(t)}{dt} =& r_{b} \cdot C(t) - r_{f} \cdot A(t) \cdot B(t) \\ \frac{dB(t)}{dt} =& r_{b} \cdot C(t) - r_{f} \cdot A(t) \cdot B(t) \\ \frac{dC(t)}{dt} =& - r_{b} \cdot C(t) + r_{f} \cdot A(t) \cdot B(t) \end{align*}\]

A simple two-compartment model could be implemented using the @reactions block via:

model = @model begin
  @param begin
    tvcl ∈ RealDomain(; lower = 0)
    tvv ∈ RealDomain(; lower = 0)
    tvka ∈ RealDomain(; lower = 0)
  end

  @pre begin
    CL = tvcl
    Vc = tvv
    Ka = tvka
  end

  @reactions begin
    Ka, Depot --> Central
    CL / Vc, Central --> 0
  end
end

The resulting dynamics can be latexified and rendered as described in Model Representation:

using Latexify
latexify(model, :reactions)

\begin{align*} \frac{dDepot(t)}{dt} =& - Ka \cdot Depot(t) \\ \frac{dCentral(t)}{dt} =& Ka \cdot Depot(t) - \frac{CL \cdot Central(t)}{Vc} \end{align*}

latexify(model, :reactions; arrow = true)

\begin{align*} \ce{ Depot &->[Ka] Central}\\ \ce{ Central &->[CL \cdot \mathrm{inv}\left( Vc \right)] \varnothing} \end{align*}

The syntax for the arrow notation could be described as rates, reactants, arrow, and products:

• arrow determines the type of reaction modeled. The most commonly used arrows are --> and <--> for uni- and bi-directional mass action reactions, respectively.
• reactants and products are essentially what you had before and what you will have after the reaction. Zero denotes the empty set and allows you to create products out of nothing (0 --> A) or to degrade reactants into nothing (A --> 0).
• rates describe the rates at which an interaction occurs. If the same expression represents more than one actual reaction (like the association and dissociation in the above example), then rates should be a tuple with one rate for each such reaction. The rates are often parameters, as defined in the @pre block, but could also be more complicated expressions involving both parameters and variables. For example, A could activate the production of B via a Michaelis-Menten function through v * A / (k + A), 0 --> B.

The parsing and interpretation of this syntax is done by Catalyst and the documentation of their DSL contains further information on how you can specify your reactions.

The main benefit of the @reactions syntax is that you can let one interaction be described in one line of code. In the common system-of-equations representation, you would often have that a single reaction gives rise to terms in multiple different equations. When you wish to modify such a reaction, you must then be careful to modify all the terms related to this reaction correctly lest your model might still run but give you an incorrect result. This mental bookkeeping that we have to do is not very burdensome in the small demos we have here but can pose real problems in larger models. The @reactions syntax may then help to keep your models nice and tidy.

@derived: Statistical modeling of observed variables in a PumasModel

This block is used to specify the assumed distributions of observed variables that are derived from the blocks above in a PumasModel. All variables are referred to as the subject's observation times which means they are vectors. This means we have to use "dot calls" on functions of dynamic variables, parameters, variables from @pre, etc.

@model begin
  @param begin
    θCL ∈ RealDomain(; lower = 0, upper = 10)
    θVc ∈ RealDomain(; lower = 0, upper = 10)
    ωη ∈ RealDomain(; lower = 0, upper = 20)
  end

  @pre begin
    CL = θCL
    Vc = θVc
  end

  @dynamics begin
    Central' = -CL / Vc * Central
  end

  @derived begin
    cp := @. Central / Vc
    dv ~ @. Normal(cp, ωη)
  end
end

We define cp (concentration in plasma) using := which means that the variable cp will not be stored in the output you get when evaluating the model's @derived block. In many cases it is easier to simply write it out like this:

@derived begin
  dv ~ @. Normal(Central / Vc, ωη)
end

This will be slightly faster. However, sometimes it might be helpful to use := for intermediary calculations in complicated expressions. An example is the proportional error model:

@model begin
  @param begin
    θCL ∈ RealDomain(; lower = 0, upper = 10)
    θVc ∈ RealDomain(; lower = 0, upper = 10)
    ωη_add ∈ RealDomain(; lower = 0, upper = 20)
    ωη_prop ∈ RealDomain(; lower = 0, upper = 20)
  end

  @pre begin
    CL = θCL
    Vc = θVc
  end

  @dynamics begin
    Central' = -CL / Vc * Central
  end

  @derived begin
    cp := @. Central / Vc
    dv ~ @. Normal(cp, sqrt(ωη_add^2 + (cp * ωη_prop)^2))
  end
end

Where we take advantage of the := line to only calculate the concentration once.

@observed: Sampled observations

If you wish to store some information from the model solution or calculate a variable based on the model solutions and parameters that has nothing to do with the statistical modeling it is useful to define these variables in the @observed block. A simple example could be that you want to store a scaled plasma concentration. This could be written like the following:

@model begin
  @param begin
    θCL ∈ RealDomain(; lower = 0, upper = 10)
    θVc ∈ RealDomain(; lower = 0, upper = 10)
    ω²η ∈ RealDomain(; lower = 0, upper = 20)
  end

  @pre begin
    CL = θCL
    Vc = θVc
  end

  @dynamics begin
    Central' = -CL / Vc * Central
  end

  @observed begin
    cp1000 = @. 1_000 * Central / Vc
  end
end

which will cause functions like simobs to store the simulated plasma concentration multiplied by a thousand (1_000).

This block is PumasModel-only.

@metadata and Descriptions

Additional details related to a PumasModel can be included in an optional @metadata block. These take the form of key/value pairs defined with = as follows where that values can take on <:Any value.

@model begin
  @metadata begin
    desc = "A short description of the model."
    timeu = u"hr"
  end
end

In the above example a textual description of the model is provided by the desc key. It must be written in plaintext and should preferably not take up more than a single line. The timeu key is assigned a Unitful value to describe the units used in the model. At present desc and timeu are the only metadata keys that are used.

In addition to the @metadata block, individual parameters, covariates, and any other named components within your models can include a short "docstring" to help readers better understand the purpose of each name. They can be added to your models in the same manner as normal docstrings used to document functions and types in normal Julia code. For example:

@model begin
  @param begin
    "Clerance (L/hr)"
    tvcl ∈ RealDomain(; lower = 0, init = 3.2)
    "Volume (L)"
    tvv  ∈ RealDomain(; lower = 0, init = 16.4)
    "Absorption rate constant (h-1)"
    tvka ∈ RealDomain(; lower = 0, init = 3.8)
    "Bioavailability"
    tvbio ∈ RealDomain(; lower = 0, init = 0.7)
    """
    Proportional RUV
    """
    σ_p ∈ RealDomain(; lower = 0, init = 0.2)
  end
  @random begin
    η ~ MvNormal(Ω)
  end
  @covariates begin
    "0  PO, 1 IV"
    formulation
    "Body Weight (kg)"
    wt
    "Age (years)"
    age
    "Race"
    racen
    "Gender"
    gender
  end
  @pre begin
    CL = tvcl * exp(η[1])
    Vc = tvv * exp(η[2])
    Ka = tvka * exp(η[3])
  end
  @dosecontrol begin
    bioav = (; Depot = tvbio * exp(η[4]))
  end
  @dynamics Depots1Central1
  @derived begin
    cp := @. Central / Vc
    "CTMx Concentrations (ng/mL)"
    dv ~ @. Normal(cp, cp*σ_p)
  end
end

For any non-RealDomain parameters you must provide a description for each element in the form of a list of descriptions, such as used to describe Ω ∈ PDiagDomain in the example above.

Always keep descriptions reasonably short, but still clear. They are used within plot labels, which have limited space available, to provide additional details.

The @emmodel macro interface

The PumasEMModel requires at least two blocks: an @error block describing the error model, and either a @param or @random block. An example minimal PumasEMModel is as follows:

@emmodel begin
  @random begin
    p ~ 1 | LogitNormal
  end
  @error begin
    y ~ Bernoulli(p)
  end
end

The possible blocks are:

• @param and @random: Defines fixed and random effects, e.g.

@random begin
  CL ~ 1 + wt | LogNormal
  θbioav ~ 1 | LogitNormal
end

Distributions specify the parameter's support. Only Normal (-∞,∞) LogNormal (0,∞), and LogitNormal (0,1) are currently supported.

These define 1 unconstrained parameter per term in the formula, for example CL = exp(μ + wt * β + η) where μ, β ∈ (-∞,∞) and wt is a covariate constant with respect to time. η is a random effect. The @param block is equivalent, with the exception that there is no random effect η. These variables are accessible in @pre, @dosecontrol, @dynamics, @post, and @error blocks:

• @param: defines population parameters without a random component.
• @random: defines population parameters with a random component.
• @covariance: defines size of the diagonal blocks of the covariance of the random effects.
• @covariates: covariates available in the @pre, @dosecontrol, @dynamics, and @post blocks.
• @pre: block evaluated before the @dynamics.
• @dosecontrol: block specifying dose control parameters. Function of @param and @random variables only.
• @init: return initial conditions for dynamical variables.
• @dynamics: dynamics, equivalent to the functionality in the @model macro.
• @post: block evaluated after the @dynamics and before the @error.
• @error: block describing the error model. Dispersion parameters are implicit. Y ~ ProportionalNormal(μ) indicates that Y has a Normal distribution with mean μ. Y must be observed subject data, while μ can can be defined in the @param, @random, @pre, @dynamics, or @post blocks.

@param: Population parameters

The @param block of a PumasEMModel defined with @emmodel has a different syntax from the @param block of a PumasModel defined with @model. For a PumasEMModel, domains are specified by a distribution with matching support: Normal for the reals, LogNormal for positive reals, and LogitNormal for those restricted to the unit interval. Only scalar distributions are supported. See the model components' documentation for more details. For example, to specify θ ∈ (0,17):

@emmodel begin
  @param begin
    unitθ ~ 1 | LogitNormal
  end
  @pre begin
    θ = 17 * unitθ
  end
  @error begin # `@error` block is mandatory for PumasEMModels
    y ~ Normal(θ)
  end
end

Unlike in a PumasModel, variables defined in the @param block are available in the @dynamics, @post, and @error blocks.

@random: Random effects

In a PumasEMModel, the @random block uses the same syntax as the @param block, and is also restricted to only the Normal, LogNormal, and LogitNormal distributions. The covariance matrix of the random effects, Ω, are defined implicitly. However, their block diagonal covariance structure can be controlled by a @covariance block. The formula definitions of the @random block also define a fixed effect per term of the formula, thus one should always define the random effects of a PumasEMModel directly in terms of the parameters that exhibit variation, rather than of the 0-mean ηs used for calculating parameters in the @pre block in a traditional PumasModel:

@emmodel begin
  @random begin
    CL ~ 1 + logwt | LogNormal   # CL = exp(log(μ_CL) + μ_CL_logwt * logwt + η_CL)
    Vc ~ 1 | LogNormal   # Vc = exp(log(μ_Vc) + η_Vc)
    F ~ 1 | LogitNormal # F  = 1/(1 + exp(-(logit(μ_F) + η_F)))
    lm ~ 1 | Normal      # lm = μ_lm + η_lm
  end
  @covariance 3, 1
  @error begin # PumasEMModels require `@error` blocks
    y ~ ProportionalNormal(CL)
  end
end

The @covariance block specifies that the covariance matrix of the four random effects (the η parameters) has a 3×3 diagonal block, shared by η_CL, η_Vc, and η_F, and a 1×1 block for η_lm. For more discussion, see the PumasEMModel-Domains documentation. Covariates appearing in the @param and @random formulas of a PumasEMModel, e.g. logwt above, do not need to be specified in the @covariates block. If one of these covariates also appears in another model block, then they must still be specified in the @covariates block, even if they appear in a formula. Covariates in a PumasEMModel formula must be constant with respect to time.

@covariance

Specify the size of the diagonal blocks in the covariance matrix of the random effects. For example:

@covariance 2 1 2

Specifies a block-diagonal covariance matrix with 2×2, 1×1 and 2×2 diagonal blocks. If unspecified, the default is to be fully diagonal, e.g. if there are 5 random effects total, the default is:

@covariance 1 1 1 1 1

@covariates

Syntax and behavior is the same as in a PumasModel. Note that covariates used inside @param and @random do not need to be specified in @covariates.

@pre

Same behavior and syntax as in a PumasModel. The @pre block is executed before the @dynamics, and the defined variables are available in each.

@init

Initial conditions for the dynamical variables. Same behavior and syntax as in a PumasModel.

@dynamics

Specify the dynamical model. Same behavior and syntax as in a PumasModel.

@post

Same behavior and syntax as the @pre block, but it is executed after the @dynamics. Because the @error block does not allow evaluating expressions or defining new variables, any variables needed in the @error block and dependent on the @dynamics should be defined in the @post block.

@error: Statistical modeling of observed variables in a PumasEMModel

This is the PumasEMModel equivalent of the @derived block in a PumasModel. The primary differences are that

1. While the @derived block applies per subject to vectors of observations across time, the @error block is mapped across time points, like the @pre and @post blocks.
2. The dispersion parameters of the @error block are defined implicitly by the specified error model.
3. The @error model does not allow for any post-processing; this must be done in the @post block.

@emmodel begin
  @random begin
    CL ~ 1 | LogNormal
    Vc ~ 1 | LogNormal
  end
  @dynamics Central1
  @post begin
    cp = Central / Vc
  end
  @error begin
    dv ~ ProportionalNormal(cp) # dv ~ Normal(cp, abs(cp)*σ)
  end
end

See the page on PumasEMModel-Error models for more information and a table of supported error models.