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 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:
using Pumas
@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:
Pumas.@model
— Macro@model
Defines a Pumas model, with the following possible blocks:
@param
: defines the model's fixed effects and corresponding domains, e.g.tvcl ∈ RealDomain(lower = 0)
.@random
: defines the model's random effects and corresponding distribution, e.g.η ~ MvNormal(Ω)
.@covariates
: Names of covariates used in the model.@pre
: pre-processes variables for the dynamic system and statistical specification.@dosecontrol
: defines dose control parameters as a function of fixed and random effects.@vars
: define variables usable in other blocks. It is recommended to define them in@pre
instead if not a function of dynamic variables.@init
: defines initial conditions of the dynamic system.@dynamics
: defines the dynamic system.@derived
: defines the statistical model of the dependent variables.@observed
: lists model information to be stored in the solution.@options
specifies model options, includingchecklinear
,inplace
, andsubject_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.
Pumas.@param
— Macro@param
Defines the model's fixed effects and corresponding domains, e.g. tvcl ∈ RealDomain(lower = 0)
. Must be used in an @model
block. For example:
@model begin
@param begin
tvcl ∈ RealDomain(lower = 0)
tvv ∈ RealDomain(lower = 0)
Ω ∈ PDiagDomain(2)
σ_prop ∈ RealDomain(lower = 0, init = 0.04)
end
end
For PumasEMModel
s, the @param
block uses a formula syntax for covariate relationships. A distribution from the normal family is used to indicate domain (currently only Normal
, LogNormal
, and LogitNormal
are supported):
@emmodel begin
@param begin
CL ~ 1 + wt | LogNormal
Vc ~ 1 | Normal
end
end
this gives CL = exp(log(θ.CL[1]) + θ.CL[2] * wt)
and VC = θ.VC
, where θ
is the parameter named tuple. If we only have 1
in the formula, then it is a scalar. If we have covariates, we have a tuple of parameters for variable.
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:
@param begin
θ ∈ RealDomain(; lower = 0, upper = 17)
end
which creates a model with a parameter that has a lower and upper bound on the allowed values.
Pumas
does not expect specific names for parameters, dependent variables, and so on. This means that fixed effects do not have to be called θ
, random effects don't have to be called η
, variability (variance-covariance) matrices for random effects don't have to be called Ω
, and so on. Pick whatever is natural for your context.
Different domains are available for different purposes. Their names and purposes are:
RealDomain
for scalar parametersVectorDomain
for vectorsPDiagDomain
for positive definite matrices with diagonal structurePSDDomain
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:
@param begin
θCL ∈ RealDomain(; lower = 0, upper = 50)
θV ∈ RealDomain(; lower = 0, upper = 500)
ω²η ∈ RealDomain(; lower = 0, upper = 20)
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:
@param begin
θ ∈ VectorDomain(2; lower = [0, 0], upper = [50, 500])
Ωη ∈ PDiagDomain(1) # no lower or upper keywords!
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 we 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 fit
ing, 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 others. 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. Note that for both of these domains, the diagonal elements are variances; and the off-diagonal elements covariances.
@param begin
θ ∈ VectorDomain(2; lower = [0, 0], upper = [50, 500])
Ωη ∈ PSDDomain(1) # no lower upper!
end
As stated above, in addition to the non-random domains, we may specify parameters via a distribution. This allows us to specify parameters 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]
# Observe that there are no commas (,) when constructing a matrix in Julia
Σ_prior = [
1.0 0.1
0.1 3.0
]
@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.
Many of the NLME model definition portions require the specification of probability distributions. The distributions in Pumas are generally defined by the Distributions
library. All of the Distributions
Distribution
types can be used throughout the Pumas model definitions. Multivariate domains defines values which are vectors while univariate domains define values which are scalars. For the full documentation of the Distribution
types, please see the Distributions
documentation.
@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.
Pumas.@random
— Macro@random
defines the model's random effects and corresponding distribution, e.g. η ~ MvNormal(Ω)
. Must be used in an @model
block. For example:
@model begin
@random begin
η ~ MvNormal(Ω)
end
end
For PumasEMModel
s, the @random
block is equivalent to the @param
s block, except that each variable defined here has an additional random effect. For example:
@emmodel begin
@random begin
CL ~ 1 + wt | LogNormal
Vc ~ 1 | Normal
end
end
this gives CL = exp(log(θ.CL[1]) + θ.CL[2] * wt + η_CL)
and VC = θ.VC + η_VC
, where θ
is the parameter named tuple. If we only have 1
in the formula, then it is a scalar. If we have covariates, we have a tuple of parameters for variable.
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:
@param begin
ωη ∈ RealDomain(; lower = 0)
end
@random begin
η ~ Normal(0, ωη)
end
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:
@param begin
α ∈ RealDomain(; lower = 0)
β ∈ RealDomain(; lower = 0)
end
@random begin
η ~ Beta(α, β)
end
It is, of course, possible to have as many univariate (or multivariate) random effects as you want. For example:
@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
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:
@param begin
ω²η ∈ VectorDomain(3, lower = 0)
end
@random begin
η ~ MvNormal(Diagonal(ω²η))
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:
@param begin
Ωη ∈ PDiagDomain(3)
end
@random begin
η ~ MvNormal(Ωη)
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:
@param begin
Ωη ∈ PSDDomain(3)
end
@random begin
η ~ MvNormal(Ωη)
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:
@param begin
ω²η ∈ RealDomain(; lower = 0)
end
@random begin
η ~ MvNormal(4, sqrt(ω²η))
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.
Supported distributions
The following table lists all distributions available for use in the @random
block. Additionally, in Bayesian models these distributions can be used as prior distributions in the @param
block.
For more details on any distribution, type ?
in the REPL followed by the distribution name, e.g. ?Beta
.
Support | Distributions |
---|---|
(0, 1) | Beta , LogitNormal |
(0, ∞) | BetaPrime , Chi , Chisq , Erlang , Exponential , Frechet , Gamma , InverseGamma , InverseGaussian , Kolmogorov , LogNormal , NoncentralChisq , Rayleigh , Weibull |
(-∞, ∞) | Cauchy , Gumbel , Laplace , Logistic , Normal , NormalCanon , NormalInverseGaussian , PGeneralizedGaussian , TDist |
Real vectors | MvNormal |
Positive vectors | MvLogNormal |
Positive definite matrices | InverseWishart , Wishart |
Correlation matrices | LKJ , LKJCholesky |
Other | Arcsine , Biweight , Cosine , Constrained , Epanechnikov , GeneralizedExtremeValue , GeneralizedPareto , Levy , LogUniform , Pareto , product_distribution , Semicircle , SymTriangularDist , Triweight , truncated , Uniform |
Univariate distributions can be scaled and shifted: Distribution μ + σ * dist
represents the distribution of a random variable of the base distribution dist
that is first scaled by σ
and then translated by μ
.
@covariates
: covariate names
The covariates in the remaining model blocks have to be specified in the @covariates
block.
Pumas.@covariates
— Macro@covariates
Defines the model's covariates, e.g. wt
. Must be used in an @model
block. For example:
@model begin
@covariates wt
end
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:
@covariates weight age OCC
or as a block:
@covariates begin
weight
age
OCC
end
The block form is mostly useful if there are a lot of covariates. Otherwise, the one-liner is preferred.
The following names are restricted from being used as covariate names: id
, amt
, time
, evid
, cmt
, rate
, duration
, ss
, ii
, route
, and tad
. These names are used for columns in the DataFrame
output form various Pumas objects, and therefore the names would clash when adding covariate information as columns.
@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.
Pumas.@pre
— Macro@pre
Pre-processes variables for the dynamic system and statistical specification. Must be used in an @model
block. For example:
@model begin
@params begin
tvcl ∈ RealDomain(lower = 0)
tvv ∈ RealDomain(lower = 0)
ωCL ∈ RealDomain(lower = 0)
ωV ∈ RealDomain(lower = 0)
σ_prop ∈ RealDomain(lower = 0, init = 0.04)
end
@random begin
ηCL ~ Normal(ωCL)
ηV ~ Normal(ωV)
end
@covariates wt
@pre begin
CL = tvcl * (wt / 70)^0.75 * exp(ηCL)
Vc = tvv * (wt / 70) * exp(ηV)
ka = CL / Vc
Q = Vc
end
end
In the @pre
block all calculations are written as if they happen at some arbitrary point in time t
. Let us see an example:
# 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
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.
In the example above, an occasion identifier or counter OCC
was used. Whenever occasions are used for Between Occasion Variability (BOV) it is important to make sure that the keyword covariates_direction
is set to :left
. Setting the value to :left
means that time-varying covariates are evaluated as the Last Observation Carried Forward (LOCF) instead of the default value of :right
that corresponds to Next Observation Carried Backward (NOCB). :left
is the default value so this issue may only arise if the setting is set to :right
- for example for NONMEM compatibility reasons.
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:
@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
where Kempt
is the so-called bile emptying rate constant that directly depends on time.
When a variable in @pre
directly depends on t
it is not possible to use analytical solutions or matrix exponential solvers. This is because the closed form solutions are derived under the assumption of piece-wise constant dynamical parameters. If you specify a closed form solution Pumas will throw an error telling you about this illegal time-dependence, and if you specified your model using written out ordinary differential equations Pumas will use a numerical integrator instead of the solver using matrix exponentials.
@dosecontrol
: Dose Control Parameters
Pumas.@dosecontrol
— Macro@dosecontrol
Define dose control parameters as a function of fixed and random effects. Options include bioav
, duration
, lags
, and rate
. Must be used in an @model
block. For example:
@model begin
@dosecontrol begin
bioav = (Depot1 = max(0, θ[5]), Depot2 = clamp(1 - θ[5], 0.0, 1.0), Central = 1)
lags = (Depot1 = 0, Depot2 = max(0.1, θ[6]), Central = 0)
end
end
The @dosecontrol
block controls the dose control parameters (DCPs) of the model: For more information, 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.
The dose control parameters are entered as NamedTuple
s. If a DCP is just set for one-compartment to have the rest default to 1.0
it is a common mistake to write rate = (Depot=θbioav)
instead of rate = (; Depot=θbioav)
. Notice the initial ;
in the second expression. This is a way to construct a NamedTuple
of length 1 in Julia.
@vars
: Short-hand notation
Suppose we have a model with a dynamic variable Central
and a volume of distribution Vc
. You can define shorthand notation for the implied plasma concentration to be used elsewhere in the model in @vars
:
Pumas.@vars
— Macro@vars
Define variables usable in other blocks. It is recommended to define them in @pre
instead if not a function of dynamic variables. Must be used in an @model
block. For example:
@model begin
@vars begin
conc = Central / Vc
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 PumasModel
s.
@init
: Initializing the dynamic system
This block defines the initial conditions of the differential equations.
Pumas.@init
— Macro@init
Defines initial conditions of the dynamic system. Must be used in an @model
block. For example:
@model begin
@init begin
Depot1 = 0.0
Depot2 = 0.0
Central = 0.0
end
end
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:
@param begin
θ ∈ VectorDomain(3; lower = [0, 0, 1], upper = [3, 1, 4])
end
@init begin
Depot = θ[2]
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.
Pumas.@dynamics
— Macro@dynamics
Defines the dynamic system by specifying differential equations. Must be used in an @model
block.
A @model
block with a @dynamics
block must not contain a @reactions
block.
Examples
A @dynamics
block with a system of differential equations:
@model begin
@dynamics begin
Depot1' = -ka * Depot1
Depot2' = ka * Depot1 - ka * Depot2
Central' = ka * Depot2 - CL / Vc * Central
end
end
A @dynamics
block with a differential equation model:
@model begin
@dynamics Depots1Central1
end
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:
@param begin
θCL ∈ RealDomain(; lower = 0, upper = 10)
θVc ∈ RealDomain(; lower = 0, upper = 10)
end
@pre begin
CL = θCL
Vc = θVc
end
@dynamics Central1
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:
@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
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 one-compartment model with first-order absorption 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
PumasModel
Parameters: tvcl, tvv, tvka
Random effects:
Covariates:
Dynamical system variables: Depot, Central
Dynamical system type: Matrix exponential
Derived:
Observed:
The resulting dynamics can be latexified and rendered as described in Model Representation:
using Latexify
latexify(model, :reactions)
\begin{align} \frac{\mathrm{d} \cdot Depot(t)}{\mathrm{d}t} &= - Ka \cdot Depot(t) \ \frac{\mathrm{d} \cdot Central(t)}{\mathrm{d}t} &= \frac{ - CL \cdot Central(t)}{Vc} + Ka \cdot Depot(t) \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 ofB
via a Michaelis-Menten function throughv * 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
.
Pumas.@derived
— Macro@derived
Defines the error model of the dependent variables. Must be used in an @model
block. For example:
@model begin
@derived begin
dv ~ @. Normal(1000 * Central / Vc, 1000 * Central / Vc * σ)
end
end
All variables are referred to at 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. The @.
macro in Julia is used for element-wise operations. It is a "dot" syntax for vectorizing function calls, which means it applies the function to each element of the array (or arrays) individually. In the context of the @derived block in Pumas, the @. macro is used to apply the operations or functions to each element of the vectors of observations across time.
@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
PumasModel
Parameters: θCL, θVc, σ
Random effects:
Covariates:
Dynamical system variables: Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dv
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)
σ ∈ 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, abs(cp) * σ)
end
end
PumasModel
Parameters: θCL, θVc, σ
Random effects:
Covariates:
Dynamical system variables: Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dv
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.
Pumas.@observed
— Macro@observed
Lists model information to be stored in the solution. Must be used in an @model
block. For example:
@model begin
@observed begin
conc = @. 1000 * Central / Vc
end
end
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)
end
@pre begin
CL = θCL
Vc = θVc
end
@dynamics begin
Central' = -CL / Vc * Central
end
@observed begin
cp1000 = @. 1_000 * Central / Vc
end
end
PumasModel
Parameters: θCL, θVc
Random effects:
Covariates:
Dynamical system variables: Central
Dynamical system type: Matrix exponential
Derived:
Observed: cp1000
which will cause functions like simobs
to store the simulated plasma concentration multiplied by a thousand (1_000
).
This block is PumasModel
-only.
@options
: Model options
When constructing the model object from the user input some checks and transformations are performed. Some of these features can be controlled or turned off by using the @options
block.
Pumas.@options
— Macro@options
Specifies model options, including the following:
checklinear::Bool
: determines whether the solver should check if the system defined in the@dynamics
block is linear. If the system is linear, setting this option totrue
(default) enables calculating the solution through matrix exponentials. If it is set tofalse
or time (t
) appears in the@pre
block, this optimization is disabled. This option can be useful when the matrix exponential solver is not superior to general numerical integrators or for debugging purposes.inplace::Bool
: controls whether the solver uses mutable operations onArray
s (true
), which can improve performance for large systems. By default, it is set totrue
for systems containing more than 5 dynamic variables, and set tofalse
otherwise.subject_t0::Bool
: determines the starting time point from which the solver integrates the system. By default, all systems are solved from time 0 to the last time point in the subject, makingsubject_t0
false
. Ifsubject_t0
is set totrue
, the solver will begin integration from the first time point found within the subject.
Must be used in an @model
block. For example (defaults are shown below):
@model begin
@options begin
checklinear = true
inplace = length(odevars) < 5
subject_t0 = false
end
end
@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.
@metadata begin
desc = "A short description of the model."
timeu = u"hr"
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
"""
Clearance (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)
"""
Variance Components: CL, Vc, Ka, bio
"""
Ω ∈ PDiagDomain(4)
"""
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, abs(cp) * σ_p)
end
end
PumasModel
Parameters: tvcl, tvv, tvka, tvbio, Ω, σ_p
Random effects: η
Covariates: formulation, wt, age, racen, gender
Dynamical system variables: Depot, Central
Dynamical system type: Closed form
Derived: dv
Observed: dv
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. A minimal example of a PumasEMModel
is as follows:
@emmodel begin
@random begin
p ~ 1 | LogitNormal
end
@error begin
y ~ Bernoulli(p)
end
end
PumasEMModel
Parameters with random effects:
p ~ (1,) | LogitNormal
The possible blocks are:
@random begin
CL ~ 1 + wt | LogNormal
θbioav ~ 1 | LogitNormal
end
Pumas.@emmodel
— Macro@emmodel
Define a PumasEMModel
. It may have the following blocks:
@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.
@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 thatY
has aNormal
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
PumasEMModel
Parameters without random effects:
unitθ ~ (1,) | LogitNormal
Unlike in a PumasModel
, variables defined in the @param
block are available in the @dynamics
, @post
, and @error
blocks.
While more complicated transforms such as the PSDDomain
aren't supported, they could be applied manually in the @pre
block. However, because the Ω
parameters are implicit, this is better handled in the @covariance
block.
@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
PumasEMModel
Parameters with random effects:
CL ~ (1, :logwt) | LogNormal
Vc ~ (1,) | LogNormal
F ~ (1,) | LogitNormal
lm ~ (1,) | Normal
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.
Parameter corresponding to 1
in the formulas specify the baseline, and are given in the parameter's natural scale. Coefficients corresponding to covariates are unconstrained. Thus in the above example, CL = exp(log(μ_CL) + μ_CL_logwt * logwt + η_CL)
, or equivalently, CL = μ_CL * exp(μ_CL_logwt * logwt) * exp(η_CL)
.
@covariance
Specify the size of the diagonal blocks in the covariance matrix of the random effects.
Pumas.@covariance
— Macro@covariance
Allows specifying a block-diagonal covariance structure for the Ω
parameters in an @emmodel
. The default is fully diagonal.
For example, to specify that the Ω
is block diagonal with 1x1
, 2x2
, and 2x2
blocks:
@emmodel begin
@random begin
Ka ~ 1 | LogNormal
CL ~ 1 | LogNormal
Vc ~ 1 | LogNormal
Q ~ 1 | LogNormal
Vp ~ 1 | LogNormal
end
@covariance 1 2 2
end
This would mean that both CL
and Vc
may be correlated with each other, and Q
and Vp
may also be, but all other pairwise correlations are zero. @covariance 5
would result in a dense matrix.
Only diagonal is currently supported when a PumasEMModel
with LaplaceI
.
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
Pumas.@covariates
— Macro@covariates
Defines the model's covariates, e.g. wt
. Must be used in an @model
block. For example:
@model begin
@covariates wt
end
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
Pumas.@pre
— Macro@pre
Pre-processes variables for the dynamic system and statistical specification. Must be used in an @model
block. For example:
@model begin
@params begin
tvcl ∈ RealDomain(lower = 0)
tvv ∈ RealDomain(lower = 0)
ωCL ∈ RealDomain(lower = 0)
ωV ∈ RealDomain(lower = 0)
σ_prop ∈ RealDomain(lower = 0, init = 0.04)
end
@random begin
ηCL ~ Normal(ωCL)
ηV ~ Normal(ωV)
end
@covariates wt
@pre begin
CL = tvcl * (wt / 70)^0.75 * exp(ηCL)
Vc = tvv * (wt / 70) * exp(ηV)
ka = CL / Vc
Q = Vc
end
end
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
Pumas.@init
— Macro@init
Defines initial conditions of the dynamic system. Must be used in an @model
block. For example:
@model begin
@init begin
Depot1 = 0.0
Depot2 = 0.0
Central = 0.0
end
end
Initial conditions for the dynamical variables. Same behavior and syntax as in a PumasModel
.
@dynamics
Pumas.@dynamics
— Macro@dynamics
Defines the dynamic system by specifying differential equations. Must be used in an @model
block.
A @model
block with a @dynamics
block must not contain a @reactions
block.
Examples
A @dynamics
block with a system of differential equations:
@model begin
@dynamics begin
Depot1' = -ka * Depot1
Depot2' = ka * Depot1 - ka * Depot2
Central' = ka * Depot2 - CL / Vc * Central
end
end
A @dynamics
block with a differential equation model:
@model begin
@dynamics Depots1Central1
end
Specify the dynamical model. Same behavior and syntax as in a PumasModel
.
@post
Pumas.@post
— Macro@post
Allows for post-processing dynamical variables, to define any variables needed in the @error
block.
Must be used in an @emmodel
block. For example:
@emmodel begin
@random begin
CL ~ 1 | LogNormal
Vc ~ 1 | LogNormal
end
@dynamics Central1
@post begin
cp = Central / Vc
end
@error begin
dv ~ CombinedNormal(cp)
end
end
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
- 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. - The dispersion parameters of the
@error
block are defined implicitly by the specified error model. - The
@error
model does not allow for any post-processing; this must be done in the@post
block.
Pumas.@error
— Macro@error
Defines the error model of the dependent variables in a PumasEMModel
. Distributions are parameterized by the mean only, the dispersion parameter(s), if any, are left implicit. See the documentation for details on the available error models. The mean parameter(s) must be (a) symbol(s). Use an earlier block, e.g. @post
to define any necessary transformations.
Must be used in an @emmodel
block. For example:
@emmodel begin
@random begin
CL ~ 1 | LogNormal
Vc ~ 1 | LogNormal
end
@dynamics Central1
@post begin
cp = Central / Vc # cp is calculated here
end
@error begin
dv ~ CombinedNormal(cp) # cp is the mean of the error distribution
end
end
See the page on PumasEMModel-Error models for more information and a table of supported error models.