@random block setup for multivariate vs univariate

1

Sometimes we need to express etas using individual distributions, especially for simulations. Is the syntax for the random block similar between multivariate Normal vs univariate? Please see below. If not, can somebody explain the implications please.

@random begin
η ~ MvNormal(Ωη)
end

@random begin
η ~ Normal(ωη) #ωη is the variance of BSV for a particular parameter
end

2

MvNormal() has only one argument whereas Normal() has two arguments. See below. If only argument is included for Normal() then the sd is assumed to be 1.0.

julia> Normal(0,1)
Normal{Float64}(μ=0.0, σ=1.0)

julia> Normal(1)
Normal{Float64}(μ=1.0, σ=1.0)

julia> Normal(0.3)
Normal{Float64}(μ=0.3, σ=1.0)

3

There are two blocks that need changes when shifting between univariate and multivariate forms:

  1. @param
    In this block, the domain of the random effect is specified, corresponding to the specification of univariate or multivariate form in the @random block. Depending on the domain used, either the standard deviations (for univariate distributions) or variances/covariances (for multivariate distributions) are given as initial estimates. There are three domains that can be used (see documentation here)
    RealDomain
    PDiagDomain
    PSDDomain

Note that initial estimates in Pumas can be specified as a separate NamedTuple outside the model block

  1. @random
    This is the block where you specify the distribution of the random variables whose corresponding initial estimates were provided in the @param block.

Here are some examples to consider:

Ex1

@param begin
    ...
    ωCL ∈ RealDomain(lower=0.0)
    ωV  ∈ RealDomain(lower=0.0)
    ...
end

param = (ωCL=0.3, ωV=0.3) # here 0.3 is the SD of the Normal Distribution

@random begin
    ηCL ~ Normal(0.0, ωCL)
    ηV  ~ Normal(0.0, ωV)
end

is equivalent to

@param begin
    ...
    Ω ∈ PDiagDomain(2)
    ...
end

param = (Ω = Diagonal([0.09 0.09]),) # here 0.09 is the diagonal element that is the variance of the MvNormal Distribution

@random begin
    η ~ MvNormal(Ω)
end

Ex2

To allow for correlations you’ll have to use MvNormal

@param begin
    ...
    Ω ∈ PSDDomain(2)
    ...
end

# e.g.
param = (Ω = diagm([0.09, 0.09]),)

# or
param = (Ω = [0.09 0.01
              0.01 0.09])

@random begin
    η ~ MvNormal(Ω)
end

Hope this helps!

1 Like