@param begin
∈ RealDomain(lower = 0, init = 5) # Typical value of Clearance
TVCL ∈ RealDomain(lower = 0, init = 50) # Typical value of Volume
TVV ∈ PDiagDomain(2) # Variance for IIV
Ω ∈ RealDomain(lower = 0, init = 0.1) # Variance for IOV of Clearance
Ω_IOV_CL ∈ RealDomain(lower = 0, init = 0.1) # Variance for IOV of Volume
Ω_IOV_Vc ∈ RealDomain(lower = 0, init = 0.1) # Proportional error
σ_prop end
@covariates OCC
@random begin
~ MvNormal(Ω)
η ~ MvNormal(Diagonal(fill(Ω_IOV_CL, 7)))
κCL ~ MvNormal(Diagonal(fill(Ω_IOV_Vc, 7)))
κVc end
@pre begin
= TVCL * exp(η[1] + κCL[OCC])
CL = TVV * exp(η[2] + κVc[OCC])
Vc end
Understanding Correlated Interoccasion Variability (IOV) Modeling in Pumas
Patrick Kofod Mogensen
1 Modeling Correlated Interoccasion Variability (IOV) in PKPD Models using 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.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.
1.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.
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.
1.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.
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.
1.3 Complete Example
1.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 * 7) # 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 aPSDDomain
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 ofDiagonal
in theMvNormal
specification - We need to input a vector of matrices to
BlockDiagonal
so we now usefill
to create a vector of 7 elements where each element is a 2x2 matrix
1.3.2 The Design
1.3.3 The Parameters
1.3.4 Simulate
1.3.5 Plot
1.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.