PumasEMModel Domains

Similar to a PumasModel, a PumasEMModel also requires you to specify domains. However, unlike a PumasModel where this is done explicitly through use of types such as VectorDomain or PSDDomain, this is done implicitly in the PumasEMModel by specifying a probability distribution. Currently, only the Normal, LogNormal, and LogitNormal distributions are supported, corresponding to domains of (-∞,∞), (0,∞), and (0,1), respectively.

Declaring parameters: @param and @random blocks

In the PumasEMModel, both the @param and @random block allow you to declare variables with with a distribution-defined domain using the same syntax. The distinction is that @param declares fixed effects, while @random implicitly adds random effects. They are declared using a formula syntax:

@param begin
    Vc ~ 1 | LogNormal
end
@random begin
    θbioav ~ 1 | LogitNormal
    CL ~ 1 | LogNormal
end

This declares one fixed effect, Vc ∈ (0,∞) and two subject-level random effects, θbioav ∈ (0,1) and CL ∈ (0,∞).

Then in simobs(model, data, param) or fit(model, data, param, approx), the parameters param must be a NamedTuple of values where the type fits in the domain. For example, the following is a valid definition of @param for the previous example, assuming our error model contained a single dv with a CombinedNormal error model:

param = (;
    Vc = 2.0,
    θbioav = 0.75,
    CL = 0.1,
    Ω = (1.0, 1.0), # diagonal elements of covariance matrix.
    σ = ((1.0, 1.0),), # (σₐ, σₚ) of `CombinedNormal`
)

Note that unlike in a PumasModel, we do not define the Ω (variance of random effects) or σ (dispersion in the error model) parameters in the @param or @random blocks, as these are defined implicitly by the number of formulas in the @random and the error models in the @error block. simobs requires we provide values for Ω and σ, although fit will initialize Ω to the identity matrix and each σ to 1.0 if absent.

The structure of Ω may be specified using a @covariance block, which allows you to specify the size of dense blocks along the diagonal. For example, if we have 6 random effects in the model,

@covariance (2, 1, 3)

indicates a 6×6 block-diagonal covariance matrix with 2×2, 1×1, and 3×3 diagonal blocks. We could initialize this via:

using StaticArrays
param = (;
    # ...
    Ω = (@SMatrix([
            1.0 0.05
            0.05 1.0
        ]), 1.3, @SMatrix([
            2.0 0.7 0.3
            0.7 1.5 0.2
            0.3 0.2 1.8
        ]))
)

Ωs are specified as a Tuple with one element per diagonal block. If a block is 1×1, then it is given by a scalar. Otherwise, by a SMatrix. If absent, Ω defaults to being diagonal.

Optionally, ω may be specified instead, where each element of the tuple corresponds to the lower triangular Cholesky factor of the corresponding block of Ω. For example, the param field using ω corresponding to Ω example above would be:

using StaticArrays
param = (
    # ...
    ω = (
        @SMatrix([
            1.0 0.0
            0.05 0.998749217771909
        ]),
        1.140175425099138,
        @SMatrix(
            [
                1.4142135623730951 0.0 0.0
                0.4949747468305832 1.1202678251204041 0.0
                0.21213203435596423 0.08480115010871586 1.3220471871080242
            ]
        )
    )
)

Formulas with covariates

It is possible to incorporate constant-time covariates in the @param and @random formulas. For example, given a logwt covariate, we can specify:

@random begin
    CL ~ 1 + logwt | LogNormal
    Vc ~ 1 | LogNormal
end

Under the hood, this defines:

CL = exp(μ_cl + logwt * β_cl + η₁)
Vc = exp(μ_vc + η₂)

where μ_cl, β_cl, and μ_vc are fixed effects, and η₁ & η₂ are the random effects corresponding to these formulas. That is, each term appearing in a formula is associated with a fixed effect. For example, 1 and logwt for CL's formula are each multiplied by a corresponding fixed effect coefficient. Formulas in the @random block also add a random effects term (η). Each of these parameters and the fixed effect are defined on the unconstrained scale, that is they're defined on (-∞,∞), and then constrained to the appropriate domain. The only valid operator is +, while valid terms are either 1 or the names of constant time covariates.

When specifying initial values, the transform/inverse transforms are applied to parameters corresponding to 1 terms. Example valid param fields for this @random block are:

param = (;
    # ...
    CL = (4.0, 0.75), # 1 + logwt
    Vc = 70.0,       # 1
    # ...
)

Thus, given logwt values of 0 and 1 respectively, the corresponding CL values are exp(log(4)) = 4 and exp(log(4) + 0.75) ≈ 8.468.