# 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,∞)`

.

Unlike in the `@param`

block of a `PumasModel`

, the distribution in an `@param`

block of a `PumasEMModel`

is only used for inferring the domain; it does not place a prior on the associated parameter. The use of a distribution here is merely to provide consistency in specification between the `@param`

and the `@random`

blocks.

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.

If fitting with SAEM, it is recommended to use larger values for `Ω`

and `σ`

than what you estimate them to be. Doing so aids exploration of the parameter space, and enables the algorithm to escape local optima.

The structure of $Ω$ may be specified using an `@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 an `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]) )
)
```

Currently fitting `PumasEMModel`

s with `LaplaceI()`

requires structurally diagonal covariance matrices.

## 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`

.