Understanding Correlated Interoccasion Variability (IOV) Modeling in Pumas

The tutorial on Modeling Interoccasion Variability (IOV) in PKPD Models using Pumas teaches how to model random effects that vary between occasions. In that situation there is a random effect κ that has one element per occasion and for each occasion it has the same variance Ω_IOV. Instead of writing out 7 different random effect definitions they were collected in one big multivariate random effect using a Diagonal matrix to represent the structure of the variance covariance matrix.

In the case where there is not one random effect that varies across occasions but several the workflow has to be adjusted slightly. This tutorial will introduce a more general approach to using a Diagonal matrix in the univariate occasion random effect case using a BlockDiagonal that allows us to represent the repeated covariance structure between the random effects.

1 Prerequisites

Before we begin, ensure that you have:

• Familiarity with the basics of nonlinear mixed-effects modeling (NLME) in Pumas including parameter and random effect specification
• A working knowledge of IOV modeling concepts in Pumas based on the previous tutorial
• The package BlockDiagonals.jl installed. using BlockDiagonals will load the package or ask you to install it if it’s not already installed.

2 How to Share Covariance Across Occasions?

When modeling IOV, it’s essential to share the variance across occasions. If there is more than one κ (interoccasion random effect) and they are not correlated we can just define them one by one.

@param begin
    TVCL ∈ RealDomain(lower = 0, init = 5)   # Typical value of Clearance
    TVV ∈ RealDomain(lower = 0, init = 50)   # Typical value of Volume
    Ω ∈ PDiagDomain(2)                       # Variance for IIV
    Ω_IOV_CL ∈ RealDomain(lower = 0, init = 0.1) # Variance for IOV of Clearance
    Ω_IOV_Vc ∈ RealDomain(lower = 0, init = 0.1) # Variance for IOV of Volume
    σ_prop ∈ RealDomain(lower = 0, init = 0.1) # Proportional error
end

@covariates OCC

@random begin
    η ~ MvNormal(Ω)
    κCL ~ MvNormal(Diagonal(fill(Ω_IOV_CL, 7)))
    κVc ~ MvNormal(Diagonal(fill(Ω_IOV_Vc, 7)))
end

@pre begin
    CL = TVCL * exp(η[1] + κCL[OCC])
    Vc = TVV * exp(η[2] + κVc[OCC])
end

Here, CL and Vc include both interindividual variability (ηCL[1],ηVc[1]) and interoccasion variability (κCL[OCC],κVc[OCC]).

The variance-covariance matrices look as follows:

\text{Cov}(κCL) = \begin{bmatrix} \Omega_{\text{IOV}_\text{CL}} & 0 & 0 & \dots & 0 \\ 0 & \Omega_{\text{IOV}_\text{CL}} & 0 & \dots & 0 \\ 0 & 0 & \Omega_{\text{IOV}_\text{CL}} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \Omega_{\text{IOV}_\text{CL}} \end{bmatrix}

and

\text{Cov}(κVc) = \begin{bmatrix} \Omega_{\text{IOV}_\text{Vc}} & 0 & 0 & \dots & 0 \\ 0 & \Omega_{\text{IOV}_\text{Vc}} & 0 & \dots & 0 \\ 0 & 0 & \Omega_{\text{IOV}_\text{Vc}} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \Omega_{\text{IOV}_\text{Vc}} \end{bmatrix}

We could also have specified this by building up the full matrix by hand. We could have made one big matrix that was 14 elements long with alternating variances of CL, Vc, CL, and so on:

\text{Cov}(κ) = \begin{bmatrix} \Omega_{\text{IOV}_\text{CL}} & 0 & 0 & \dots & 0 \\ 0 & \Omega_{\text{IOV}_\text{Vc}} & 0 & \dots & 0 \\ 0 & 0 & \Omega_{\text{IOV}_\text{CL}} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \Omega_{\text{IOV}_\text{Vc}} \end{bmatrix}

Then the model would be

@param begin
    TVCL ∈ RealDomain(lower = 0, init = 5)   # Typical value of Clearance
    TVV ∈ RealDomain(lower = 0, init = 50)   # Typical value of Volume
    Ω ∈ PDiagDomain(2)                       # Variance for IIV
    Ω_IOV_CL ∈ RealDomain(lower = 0, init = 0.1) # Variance for IOV of Clearance
    Ω_IOV_Vc ∈ RealDomain(lower = 0, init = 0.1) # Variance for IOV of Volume
    σ_prop ∈ RealDomain(lower = 0, init = 0.1) # Proportional error
end

@covariates OCC

@random begin
    η ~ MvNormal(Ω)
    κ ~ MvNormal(Diagonal(repeat([Ω_IOV_CL, Ω_IOV_Vc], 7)))
end

@pre begin
    CL = TVCL * exp(η[1] + κ[OCC*2-1])
    Vc = TVV * exp(η[2] + κ[OCC*2])
end

This works because the nth occasion indexed by OCC will then pick the CL component at index 1,3,5,7,…,13 and the Vc component will be at 2,4,6, …, 14.

2.1 Introducing Covariances

Now that multiple interoccasion variables are introduced we can expand the framework to include covariances

\text{Cov}(κ) = \begin{bmatrix} \Omega_{\text{CL}} & \Omega_{\text{CL,Vc}} & 0 &0 & \dots & 0& 0 \\ \Omega_{\text{CL,Vc}} & \Omega_{\text{Vc}} & 0 & 0 &\dots & 0& 0 \\ 0 & 0 & \Omega_{\text{CL}} & \Omega_{\text{CL,Vc}}& \dots & 0& 0 \\ 0 & 0 & \Omega_{\text{CL,Vc}} & \Omega_{\text{Vc}} & \dots & 0& 0 \\ \vdots & \vdots & \vdots & \vdots &\ddots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & \dots & \Omega_{\text{CL}}& \Omega_{\text{CL,Vc}} \\ 0 & 0 & 0 & 0 & \dots & \Omega_{\text{CL,Vc}}& \Omega_{\text{Vc}} \end{bmatrix}

Here it is clear that the variance-covariance matrix is just a repetition of the 2x2 matrix specifying the per occasion variance-covariance structure just like the repeated variance diagonal in the example above. In Pumas such a matrix can be constructed easily using the BlockDiagonal type.

using BlockDiagonals
V = [1 0.2; 0.2 2]
BV = BlockDiagonal(fill(V, 7))

14×14 BlockDiagonal{Float64, Matrix{Float64}}:
 1.0  0.2  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.2  2.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.2  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.2  2.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.2  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.2  2.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.2  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.2  2.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.2  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.2  2.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.2  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.2  2.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.2
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.2  2.0

The variable V above is the variance-covariance matrix with variances 1 and 2 and covariance 0.2. This could represent the clearance and volume IOVs from above. The BlockDiagonal type takes one of these matrices and repeats it down the diagonal such that the final interpretation is exactly as above: CL entries of the κ are in 1,3,…,13 and the Vc entries are in 2,4,…,14.

3 Complete Example

3.1 The model

Putting it all together, here’s the full model code:

using Pumas, BlockDiagonals

IOV_PK_model = @model begin
    @param begin
        TVCL ∈ RealDomain(lower = 0, init = 5)   # Typical value of Clearance
        TVV ∈ RealDomain(lower = 0, init = 50)   # Typical value of Volume
        Ω ∈ PDiagDomain(2)                       # Variance for IIV
        Ω_IOV ∈ PSDDomain(2)                     # Variance-Covariance for IOV
        σ_prop ∈ RealDomain(lower = 0, init = 0.1) # Proportional error
    end

    @covariates OCC

    @random begin
        η ~ MvNormal(Ω)
        κ ~ MvNormal(BlockDiagonal(fill(Ω_IOV, 7)))
    end

    @pre begin
        CL = TVCL * exp(η[1] + κ[OCC*2-1])
        Vc = TVV * exp(η[2] + κ[OCC*2])
    end

    @dynamics Central1

    @derived begin
        cp = @. Central / Vc
        dv ~ @. Normal(cp, abs(cp) * σ_prop)
    end
end

PumasModel
  Parameters: TVCL, TVV, Ω, Ω_IOV, σ_prop
  Random effects: η, κ
  Covariates: OCC
  Dynamical system variables: Central
  Dynamical system type: Closed form
  Derived: cp, dv
  Observed: cp, dv

The key things to notice are:

• Ω_IOV is now a PSDDomain so the input to any simulation or fitting procedure must be a 2x2 matrix. This is true even if there are more or less occasions as we only parameterize each block we do not estimate the entire variance covariance matrix as if all elements are free.
• We use BlockDiagonal instead of Diagonal in the MvNormal specification
• We need to input a vector of matrices to BlockDiagonal so we now use fill to create a vector of 7 elements where each element is a 2x2 matrix

3.2 The Design

pop1 = map(
    subj -> Subject(
        id = subj,
        events = DosageRegimen(100, addl = 6, ii = 24),
        covariates = (; OCC = [1, 2, 3, 4, 5, 6, 7]),
        covariates_time = [0, 24, 48, 72, 96, 120, 144],
    ),
    1:10,
)

Population
  Subjects: 10
  Covariates: OCC
  Observations: 

3.3 The Parameters

params = (
    TVCL = 4,
    TVV = 10,
    Ω = Diagonal([0.04, 0.04]),
    Ω_IOV = [0.4 0.01; 0.01 0.5],
    σ_prop = 0.2,
)

(TVCL = 4,
 TVV = 10,
 Ω = [0.04 0.0; 0.0 0.04],
 Ω_IOV = [0.4 0.01; 0.01 0.5],
 σ_prop = 0.2,)

3.4 Simulate

using Random
Random.seed!(1234)
sims = simobs(IOV_PK_model, pop1, params, obstimes = 0:1:192)

Simulated population (Vector{<:Subject})
  Simulated subjects: 10
  Simulated variables: cp, dv

3.5 Plot

using PumasUtilities
sim_plot(sims, observations = [:cp])

4 Conclusion

Modeling interoccasion variability in Pumas is straightforward once you understand the concepts:

• Shared Variance-Covariance: By sharing a variance-covariance matrix Ω_IOV across occasions, we ensure that the variability and correlation is the same across occasions.
• Indexing: Using the occasion covariate (OCC) to index several IOV random effects allows each occasion to have unique effects while maintaining shared distributional assumptions.
• Implementation: Pumas provides a flexible framework to incorporate IOV into your models with minimal additional code.

Note: Always ensure that your data includes an appropriate occasion covariate and that it’s correctly specified in your Pumas model. This is especially true when introducing covariances as there is typically not a lot of data in each occasion to estimate these.