How to handle several random effects with pairwise correlation between some

How-to
1

Sometimes it is useful to have some random effects correlated and others uncorrelated with each other. In Pumas, this is possible by specifying separate random effects for those that are correlated and those that are not correlated. In the below example, we have that the clearance and volume of distribution parameters are correlated through the first and second component of the η_vc_cl random effect, while the other random effect on the rate of absorption η_ka is uncorrelated with the other two.

model = @model begin
    @param begin
        tvvc ∈ RealDomain(lower = 0)
        tvcl ∈ RealDomain(lower = 0)
        tvka ∈ RealDomain(lower = 0)
        Ωvc_cl ∈ PSDDomain(2)
        ωka ∈ RealDomain(lower = 1e-5)
        σ_pk ∈ RealDomain(lower = 1e-5)
    end
    @random begin
        η_vc_cl ~ MvNormal(Ωvc_cl)
        η_ka ~ Normal(0, ωka)
    end
    @pre begin
        Ka = tvka*exp(η_ka)
        CL = tvcl*exp(η_vc_cl[1])
        Vc = tvvc*exp(η_vc_cl[2])
    end
    @dynamics Depots1Central1
    @derived begin
        conc_pk ~ @. LogNormal(log(Central/Vc), σ_pk)
    end
end

doses = DosageRegimen(10, addl = 3, ii = 12)

raw_subject = Subject(id="1", events = doses)

param = (tvvc = 2.1, tvcl = 2.1, tvka = 0.9, Ωvc_cl = [0.4 0.2; 0.2 0.4], ωka = 0.3, σ_pk = 0.2)

subject = Subject(simobs(model, raw_subject, param; obstimes = 0:0.5:72))

using PumasUtilities
sim_plot(subject)

The simulation gives a trajectory as follows

Skærmbillede 2023-09-17 kl. 21.30.02
Skærmbillede 2023-09-17 kl. 21.30.021578×1178 77 KB

To get a confirmation that the two components are indeed correlated, we can simulate many subjects as follows

raw_population = [Subject(id=string(id), events = doses) for id = 1:1000]
sim_population = simobs(model, raw_population, param; obstimes = 0:0.5:72)
population = Subject.(sim_population)
sim_plot(population)

and this gives trajectories as follows

Skærmbillede 2023-09-17 kl. 21.32.42
Skærmbillede 2023-09-17 kl. 21.32.421578×1178 141 KB

Then, we can look at their tabular representation using the DataFrame constructor

sim_population_df = DataFrame(sim_population)
# First the first row of each subject
sim_popuplation_df_first = unique(sim_population_df, [:id])

and from this we can look at correlations between the different random effects using an AlgebraOfGraphics plot

draw(data(sim_population_df) * (visual(AlgebraOfGraphics.Scatter) + AlgebraOfGraphics.linear()*visual(color=:red) * mapping([:η_vc_cl_1, :η_vc_cl_2,:η_ka_1 ], [:η_vc_cl_1  :η_vc_cl_2 :η_ka_1 ], col = dims(1), row = dims(2)))

Skærmbillede 2023-09-17 kl. 21.48.40
Skærmbillede 2023-09-17 kl. 21.48.401578×1178 74.9 KB

and we see that the correlations are as expected.

This example had a pair with correlations and a third random effect that was uncorrelated with the other two, but the scalar random effect for rate of absorption that had a Normal distribution could just as well have been multivariate, or we could have had ever more random effects each with their own correlation structure between the elements such as the following pseudo-example

@model begin
    @param begin
        [...]
         Ω1 ∈ PSDDomain(2)
         Ω2 ∈ PSDDomain(4)
         Ω3 ∈ PSDDomain(3)
        [...]
    end
    @random begin
        η1 ~ MvNormal(Ω1)
        η2 ~ MvNormal(Ω2)
        η3 ~ MvNormal(Ω3)
    end
    [...]
end

Hope this helps!

2 Likes